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Abstract 

This  thesis  addresses  the  problem  of  optimally  selecting  and  specifying  satel¬ 
lite  payloads  for  inclusion  on  a  satellite  bus  to  be  launched  into  a  constellation. 
The  objective  is  to  select  and  specify  payloads  so  that  the  total  lifetime  utility  of 
the  constellation  is  maximized.  The  satellite  bus  is  limited  by  finite  power,  weight, 
volume,  and  cost  constraints.  This  problem  is  modeled  as  a  classical  knapsack  prob¬ 
lem  in  one  and  multiple  dimensions,  and  dynamic  programming  and  binary  integer 
programming  formulations  are  provided  to  solve  the  problem.  Due  to  the  compu¬ 
tational  complexity  of  the  problem,  the  solution  techniques  include  exact  methods 
as  well  as  four  heuristic  procedures  including  a  greedy  heuristic,  two  norm-based 
heuristics,  and  a  simulated  annealing  heuristic.  The  performance  of  the  exact  and 
heuristic  approaches  is  evaluated  on  the  basis  of  solution  quality  and  computation 
time  by  solving  a  series  of  notional  and  randomly-generated  problem  instances.  The 
numerical  results  indicate  that,  when  an  exact  solution  is  required  for  a  moderately- 
sized  constellation,  the  integer  programming  formulation  is  most  reliable  in  solving 
the  problem  to  optimality.  However,  if  the  problem  size  is  very  large,  and  near- 
optimal  solutions  are  acceptable,  then  the  simulated  annealing  algorithm  performs 
best  among  the  heuristic  procedures. 
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OPTIMIZING  MEAN  MISSION  DURATION 
FOR  MULTIPLE-PAYLOAD  SATELLITES 


1.  Introduction 

1 . 1  Background 

The  U.S.  Department  of  Defense  (DoD)  spends  approximately  $18  billion  dol¬ 
lars  annually  on  the  development,  procurement,  and  operation  of  satellites  and  other 
space-based  assets.  According  to  a  2003  General  Accounting  Office  report  [9] ,  many 
of  these  acquisition  programs  consistently  suffer  from  cost  overruns  and  delays.  In¬ 
dividual  satellites  and  satellite  constellations  constitute  a  sizable  portion  of  these 
space-based  systems.  However,  it  is  possible  that  the  costs  of  the  satellite  acquisi¬ 
tion  process  may  be  greatly  reduced  through  the  application  of  methodologies  that 
more  effectively  allocate  resources. 

Satellites  range  in  complexity  from  relatively  simple  $40-million  GPS  naviga¬ 
tion  satellites  to  complex  and  expensive  half-billion  dollar  MILSTAR  tactical  com¬ 
munication  satellites.  Although  the  current  satellite  design  trend  is  toward  simple 
and  inexpensive  satellites  clue  to  mitigation  of  launch  and  component  failures,  pre¬ 
vailing  circumstances  still  necessitate  the  design  and  use  of  expensive,  complex,  mul¬ 
tipurpose  satellites.  Because  of  the  nature  of  funding,  organizations  often  receive 
large,  irregular  monetary  allocations  to  acquire  satellite  systems,  and  in  order  to 
use  funding  more  effectively,  satellites  are  designed  to  have  as  many  mission  capa¬ 
bilities  as  possible.  Furthermore,  multiple  satellites  often  work  in  concert  forming  a 
constellation.  As  satellites  in  the  constellation  begin  to  degrade  and  lose  function¬ 
ality,  decisions  about  how  to  best  replace  them  given  limited  financial  and  material 
resources  must  be  made.  Deciding  which  mission  capabilities  (i.e.  payloads)  to  in- 
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elude  on  new  satellites  is  paramount.  Payloads  must  be  selected  according  to  the 
functions  they  perform  and  given  both  physical  and  design  parameter  specifications 
to  ensure  they  are  compatible  with,  and  operate  effectively  on,  the  satellite  bus. 

One  of  the  critical  specifications  of  a  satellite  payload  is  its  mean  mission  du¬ 
ration  (MMD).  MMD  is,  approximately,  a  measure  of  the  duration  mission  planners 
can  expect  a  payload  to  be  functional.  MMD  is  not  a  mean  in  the  mathematical 
sense,  nor  does  it  denote  how  long  a  payload  will  actually  function.  It  simply  pro¬ 
vides  mission  planners  a  reasonable  expectation  of  the  amount  of  time  the  payload 
will  be  useful.  As  shown  in  Figure  1.1,  MMD  is  usually  less  than  the  design  life  of 
the  payload.  Increasing  a  payload’s  MMD  is  analogous  to  increasing  its  reliability. 
This  is  accomplished  by  adding  redundant  systems  and  components  as  well  as  using 
materials  less  susceptible  to  degradation.  Such  measures  increase  the  payload’s  cost, 
weight,  volume,  and  other  requirements.  These  are  resources  for  which  a  satellite 
bus  has  only  a  finite  capacity. 


Survival 

Probability 


Figure  1.1  Graphical  depiction  of  payload  MMD. 


Some  methodologies  currently  exist  for  selecting  and  specifying  satellite  pay- 
loads.  However,  they  are  either  very  general  and  qualitative  or  so  specific  that  their 
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application  to  general  specification  problems  is  very  limited.  This  research  seeks  to 
develop  an  analytical  methodology  to  select  and  specify  satellite  payloads  that  has 
sufficient  generality  for  application  to  virtually  any  type  of  satellite  constellation. 
Such  a  method  would  take  quantifiable  characteristics  present  in  every  type  of  satel¬ 
lite  and  payload  and  use  that  information  in  a  methodical  approach  to  make  payload 
selection  and  specification  decisions.  The  use  of  operations  research  techniques  can 
greatly  reduce  the  cost  of  a  satellite  constellation  by  allocating  resources  more  ef¬ 
fectively  while  simultaneously  increasing  its  ability  to  achieve  the  overall  mission 
objectives. 

1.2  Problem  Definition  and  Methodology 

Consider  a  satellite  constellation  observed  at  fixed  times.  Each  satellite  has  a 
set  of  payloads  associated  with  it,  and  at  each  observation,  the  functional  status  of 
each  payload  in  the  constellation  is  known.  Also  known  are  the  MMDs  and  time  in 
service  of  payloads  initially  in  the  constellation.  The  survival  distributions  of  similar 
or  dissimilar  payloads  are  assumed  to  be  statistically  independent.  At  predetermined 
times,  single-satellite  buses  will  be  equipped  with  payloads  from  a  fixed  set  of  all 
available  payloads  and  launched  into  the  constellation.  Any  selected  payload  must 
be  assigned  a  MMD  specification  from  a  finite  set  of  MMD  specifications  available 
to  that  payload.  Each  bus  has  finite  power,  cost,  weight,  and  volume  constraints.  A 
payload’s  type  denotes  its  specific  function  and  is  independent  of  its  MMD  specifica¬ 
tion.  Moreover,  each  payload  has  a  utility  associated  with  it  though  utility,  as  it  is 
used  in  this  thesis,  does  not  satisfy  the  strict  definition  of  utility.  Payload  utility  is 
assumed  to  be  a  function  of  the  payload’s  relative  importance,  MMD  specification, 
and  the  expected  number  of  functional  like  payloads  in  the  constellation.  The  impor¬ 
tance  of  a  payload  is  analogous  to  its  value  to  the  mission.  It  is  only  dependent  on 
the  payload’s  type  and  is  not  affected  by  the  payload’s  MMD  specification.  Utility 
dependence  on  the  number  of  like-type,  functioning  payloads  allows  a  payload’s  util- 
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ity  function  to  model  diminishing  marginal  returns.  It  is  assumed  that  a  payload’s 
power  consumption  is  proportional  to  its  utility.  For  a  given  payload  type  and  MMD 
specification,  discrete  functions  are  assumed  to  exist  that  give  the  cost,  weight,  and 
volume  resources  consumed  by  the  payload. 

The  objective  is  to  select  and  specify  payloads  for  the  satellites  being  launched 
such  that  the  total  overall  utility  of  the  constellation  is  maximized.  Payload  utility 
can  be  static  (constant  over  time)  or  dynamic  (changing  over  time).  Additionally,  it 
can  be  deterministic  (known  with  certainty)  or  stochastic  (probabilistic).  Enumer¬ 
ating  each  combination,  utility  can  be  characterized  as  one  of  the  following:  static 
and  deterministic,  dynamic  and  deterministic,  static  and  stochastic,  or  dynamic  and 
stochastic.  The  characterization  of  utility  is  integral  to  the  form  of  the  resulting 
mathematical  model  for  payload  selection.  The  problem  of  selecting  and  specifying 
satellite  payloads  is  similar  to  a  class  of  mathematical  programming  problems  known 
as  knapsack  problems.  Given  a  set  of  items,  each  having  an  associated  profit  and 
weight,  the  knapsack  problem  seeks  to  place  them  in  a  knapsack  of  finite  weight- 
capacity  such  that  the  profit  of  included  items  is  maximized.  A  payload  selection 
model  for  a  single  satellite  will  be  developed  for  each  of  the  four  characterizations 
of  payload  utility  and  shown  to  be  a  relaxation  of  a  type  of  multidimensional  knap¬ 
sack  problem.  The  dynamic  and  stochastic  case  most  realistically  describes  the 
nature  of  actual  payload  utility  and  is  extended  to  a  multi-satellite  case.  A  so¬ 
lution  methodology  will  be  developed  to  solve  the  multi-satellite  problem.  Exact 
solutions  to  knapsack  problems  generally  require  either  dynamic  programming  or 
integer  programming  formulations.  Therefore,  in  order  to  apply  knapsack-based  so¬ 
lution  techniques  to  the  multi-satellite  payload  selection  problem,  the  multi-satellite 
problem  will  be  formulated  as  both  a  dynamic  program  and  an  integer  program. 

A  dynamic  program  can  be  solved  (at  great  expense)  by  completely  enumerat¬ 
ing  the  state  space,  and  such  a  method  can  be  applied  to  the  dynamic  programming 
formulation  of  the  payload  selection  problem.  Integer  programs  are  solved  primarily 
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using  branch-and-bound  or  branch- and-cut  algorithms.  Commercial  IP  solvers,  like 
Xpress,  apply  these  algorithms  in  concert  with  preprocessing  and  search  heuristics  to 
attain  solutions  more  quickly.  To  provide  a  solution  to  the  integer  programming  for¬ 
mulation,  the  Xpress  solver  will  be  used.  Moreover,  heuristic  solution  methods  will 
be  developed.  An  elementary  heuristic  for  the  one-dimensional  knapsack  problem 
is  based  on  the  profit-to-cost  ratio  of  the  items.  This  can  be  extended  to  solve  the 
relaxed,  multidimensional  knapsack  problem  associated  with  the  payload  selection 
problem.  Furthermore,  the  ability  of  a  more  advanced  heuristic  like  simulated  an¬ 
nealing  to  solve  the  payload  selection  problem  will  also  be  explored.  Exact  methods 
guarantee  the  optimality  of  the  solution  but  can  take  considerable  time.  Heuristics 
are  generally  faster;  however,  they  provide  no  guarantee  of  convergence  to  optimal¬ 
ity.  It  is  desirable  to  determine  which  method  can  best  be  applied  to  the  problem 
of  satellite  payload  selection  and  specification.  The  performance  of  the  exact  and 
heuristic  solution  methods  will  be  compared  by  applying  them  to  multiple  problem 
instances  using  notional  and  randomly-generated  payload  and  satellite  data  sets. 
Ultimately,  one  of  the  main  objectives  of  this  thesis  is  to  provide  recommendations 
for  which  technique  performs  best  to  provide  optimal,  or  near-optimal,  solutions  in 
a  reasonable  amount  of  time. 

1.3  Thesis  Outline 

The  remainder  of  the  thesis  is  organized  as  follows.  Chapter  2  presents  a  review 
of  literature  associated  with  payload  selection  methodologies  as  well  as  descriptions 
and  main  results  for  knapsack  problems.  In  chapter  3,  the  assumptions  and  math¬ 
ematical  models  for  both  the  single-satellite  and  multi-satellite  payload  selection 
problems  are  discussed  followed  by  dynamic  programming  and  integer  programming 
formulations  of  the  multi-satellite  case.  Chapter  4  introduces  exact  and  heuristic 
solution  methods  for  solving  the  payload  selection  problem  as  well  as  prospective 
functional  forms  for  the  payload  utility  and  survival  functions.  The  solution  methods 
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are  subsequently  applied  to  problem  instances  using  both  notional  and  randomly- 
generated  data  sets  for  payload  resource  requirements  and  satellite  resource  capac¬ 
ities.  Finally,  chapter  5  summarizes  conclusions,  unresolved  issues,  and  possibilities 
for  future  extensions  of  the  research. 
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2.  Relevant  Literature 


This  chapter  reviews  the  literature  pertaining  to  the  selection  and  MMD 
specification  of  satellite  payloads.  The  chapter  begins  with  a  review  of  current  pay- 
load  design  and  selection  methodologies  and  proceeds  to  introduce  a  closely-related 
problem  known  as  the  knapsack  problem.  Because  the  knapsack  problem  provides  a 
framework  for  payload  selection,  potential  exact  and  heuristic  solution  methods  for 
the  knapsack  problem  are  also  discussed. 

2.1  Payload  Selection  Strategies 

The  literature  on  payload  selection  and  specification  is  relatively  sparse.  Most 
literature  concerning  satellite  payloads  addresses  the  intricacies  associated  with  their 
design  and  construction.  Larson  and  Wertz  [13]  present  a  ten-step  methodology  for 
payload  selection  and  specification  based  on  a  satellite’s  mission  objectives.  The  pro¬ 
cess  begins  by  translating  the  relatively  general  mission  objectives  of  a  satellite  into 
payload  objectives  that  denote  the  specific  functions  payloads  are  required  to  per¬ 
form.  Next,  subject  trades  are  conducted  in  which  subjects ,  the  specific  objects  with 
which  payloads  interact,  are  identified.  Together  the  subjects  and  payload  objectives 
allow  generation  of  a  payload  operations  concept.  The  payload  operations  concept 
identifies  what  is  necessary  to  enable  payloads  to  both  perform  their  functions  and 
communicate  results. 

Next,  throughput  and  performance  requirements  of  the  payloads  can  be  iden¬ 
tified  followed  by  the  actual  identification  of  candidate  payloads.  Because  tasks  can 
often  be  broken  down  or  shared  in  many  different  ways,  many  possibilities  exist  for 
candidate  payloads.  The  characteristics  of  each  candidate  payload  must  be  estimated 
including  its  performance  and  resource  requirements.  Once  the  characteristics  are 
estimated,  a  base-line  of  prospective  payloads  is  selected.  This  selection  is  usually 
based  on  cost  versus  performance  considerations.  Next,  both  life-cycle  cost  and  op- 
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erability  need  to  be  assessed.  Although  this  step  ultimately  seeks  to  define  payload 
utility  as  a  function  of  cost,  it  is  very  complex  and  involves  extensive  contact  with 
the  satellite’s  prospective  users  and  possible  mission  objective  tradeoffs.  Finally,  the 
requirements  imposed  by  the  payloads  on  the  satellite  bus,  ground  operations,  and 
mission  control  are  enumerated.  The  process  is  concluded  by  thoroughly  document¬ 
ing  all  conclusions  and  may  be  repeated  multiple  times. 

Although  this  selection  methodology  provides  a  general  framework  for  pay- 
load  selection  and  specification,  it  does  not  provide  a  specific  quantitative  method 
for  decisions  regarding  the  inclusion,  exclusion,  or  specification  of  payloads.  Bell 
[2]  presents  a  methodology  applied  to  secondary  payloads  on  GPS  satellites.  A 
secondary  nuclear  detonation  detection  system  (NDS)  payload  is  included  on  ev¬ 
ery  GPS  satellite.  However,  NDS  ground  systems  can  only  monitor  twenty-four  of, 
at  that  time,  twenty-nine  GPS  satellites.  Determining  which  satellites  to  monitor 
is  equivalent  to  a  payload  selection  problem  in  which  payloads  are  to  be  selected 
for  twenty-four  satellites.  The  decision  of  which  satellites  to  monitor  was  based  on 
their  individual  value.  A  satellite  is  valued  by  its  contribution  to  coverage ;  that  is, 
its  ability  to  both  observe  and  detect  a  nuclear  detonation.  A  satellite’s  real-time 
connectivity,  optical  sensor  sensitivity,  and  orbital  location  were  used  to  compute  a 
contribution  to  coverage  coefficient.  The  first  step  of  the  solution  involves  determin¬ 
ing  an  initial  coverage  that  ensures  a  relatively  even  distribution  of  satellites  over 
each  orbital  plane.  A  heuristic  is  applied  iteratively  that  replaces  spare  satellites  in 
each  plane.  For  each  successive  solution,  the  coverage  of  the  entire  constellation  is 
calculated  using  a  classified,  GPS/MS  simulation  program. 

2.2  Potential  Optimal  Solution  Methodologies 

The  satellite  payload  selection  problem  seeks  to  maximize  the  total  overall 
utility,  or  some  other  benefit,  subject  to  resource  constraints.  This  is  analogous  to 
a  class  of  problems  known  as  knapsack  problems.  This  section  discusses  both  the 
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one-  and  multidimensional  knapsack  problems  (KP/MKP)  as  well  as  exact  solution 
methods  for  solving  them. 

2.2.1  Overview  of  Knapsack  Problems 

The  one-dimensional  knapsack  problem  is  one  of  the  oldest  and  most  widely- 
studied  mathematical  programming  problems.  Given  a  knapsack  of  finite  capacity  b, 
which  items,  each  having  weight  al  and  profit  Ci,  i  =  1,2, ...  ,1,  should  be  placed  in 
the  knapsack  such  that  total  benefit  is  maximized  and  the  capacity  of  the  knapsack 
is  not  exceeded?  Let  xt  —  1  denote  the  inclusion  of  item  i  and  Xi  =  0  denote  the 
exclusion  of  item  i,  i  —  1,2, ...  ,1.  Stated  mathematically  it  is: 

i 

Maximize  E  CiXi 

i=  1 
l 

Subject  to  aiXi  <  b 
1=1 

Xi  g  {0,1} 

di,  b  G  Z+,  i  —  l,2,...,l. 

The  knapsack  problem  is  NP-hard,  and  no  algorithm  is  known  that  yields 
an  optimal  solution  in  polynomial-time;  however,  fully  polynomial  approximation 
algorithms  do  exist  [17].  Many  variations  of  the  one- dimensional  KP  have  been 
formulated  to  include  the  multiple-choice  KP  in  which  the  item  set  0  is  partitioned 
into  k  subsets,  @i,  02, . . . ,  0*,,  such  that  only  one  item  in  each  subset  may  be  selected. 
Another  variation  is  the  bounded  KP  in  which  each  Xi  G  Z+  and  is  bounded  by 
bi  G  Z+  such  that  Xi  <  b, .  Despite  their  differences,  both  the  multiple-choice  and 
bounded  KPs  are  NP-hard. 

Many  important  results  of  the  one-dimensional  KP  center  around  its  LP- 
relaxation.  Dantzig’s  [4]  classical  method  to  solve  the  LP-relaxation  is  to  order 
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objects  by  decreasing  profit-to-weight  ratio,  where 


Cl 

C2 

Cl 

>  —  >  . 

U\ 

a  2 

a,[ 

(2.1) 


The  solution  proceeds  with  the  insertion  of  item  1,  and  continues  by  inserting 
items  into  the  knapsack  by  successively  decreasing  profit-to-weight  ratios.  Eventu¬ 
ally,  item  s  will  be  reached  which  cannot  fit  into  the  knapsack,  ft  can  be  shown 
[17]  that  the  optimal  solution  to  the  LP-relaxation  is  given  by  xs  =  b/as,  where 

_  5—  1 

b  =  b  —  ^2  a.i.  Item  s  is  deemed  the  break  or  critical  item.  Simply  including  the  first 

i=  1 

s  —  1  items  in  the  knapsack  often  provides  a  good  solution  and  forms  the  basis  of  a 
greedy,  profit-to-cost  ratio  heuristic.  Based  on  integrality  of  all  profits  and  weights, 
Dantzig  [4]  showed  an  upper-bound  of  the  knapsack  solution  is 


S—  1 

+ 

i= 1 


(2.2) 


This  upper  bound  plays  a  role  in  branch-and-bound  algorithms  discussed  later  in 
this  chapter. 

Consider  a  generalization  of  the  knapsack  problem  to  the  multidimensional 
case.  The  multidimensional  knapsack  problem  assumes  that  the  knapsack  has  mul¬ 
tiple  constraints  akin  to  multiple  dimensions.  For  a  problem  with  m  constraints,  the 
mathematical  program  for  the  MKP  is: 

i 

Maximize  E  CjXj 

j= i 

i 

Subject  to  aijXj  <  bt ,  i  =  1,2,...,  m 

3= 1 

Xj  e  {0,1},  j  =  1,2, ...,/, 
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where  at]  is  the  size  of  object  j  on  the  ith  dimension.  The  MKP  is  NP-hard;  however, 
unlike  the  one-dimensional  KP,  no  fully-polynomial  approximation  algorithms  exist 
unless  P  =  NP  [17];  however,  polynomial  approximations  do  exist  [8].  Variations  for 
the  MKP  also  exist  including  the  multi-choice,  multidimensional  knapsack  problem 
(MCMKP).  A  formulation  for  the  MCMKP  is  provided  for  this  thesis  based  on  an 
adaptation  of  Martello  and  Toth’s  [17]  formulation  for  the  one- dimensional,  multi¬ 
choice  knapsack  problem.  Define  the  following  quantities: 


0 

01,  ©2,  •  •  •  ,  ©fc 

Cj 

Ciij 

h 


Item  set 

Partition  of  item  set  0 

k 

where,  0*  =  0  and  0,  fl  0^  =  0,  i  7^  j 
1=1 

Benefit  of  item  j,  j  —  1,2, ...  ,1 

Size  of  item  j  on  dimension  % 

j  =  1,2,...,/,  i  =  1,2,  ...,D 

Limit  of  knapsack  dimension  i,  i  =  1,2, . . . ,  D 


Let 


{1,  if  item  j  is  included 
0,  if  item  j  is  not  included 


The  MCMKP  may  be  formulated  as 


1 


Maximize 

cixi 
3= 1 

(2,3) 

Subject  to 

1 

aijxj 

7=1 

<  bh  7  =  1,2,  ...,D 

(2,4) 

=  1,  i  —  1,2, . . .  ,k 

(2,5) 

Xj 

e  {0,1},  j  —  1,2, ...  ,l. 

(2,6) 
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Many  different  algorithms  have  been  developed  to  solve  the  one-dimensional 
and,  to  a  lesser  extent,  the  multidimensional  knapsack  problems  exactly.  These 
methods  are  primarily  based  on  dynamic  or  integer  programming  formulations.  Each 
is  addressed  in  turn  in  the  following  subsections. 

2.2.2  Dynamic  Programming  and  KP/MKP 

Dynamic  programming  provides  a  method  to  analyze  sequential  decision  pro¬ 
cesses  [5].  It  is  based  on  a  recursive  relationship  known  as  Bellman’s  functional 
equation.  In  general,  dynamic  programming  requires  a  state  space  denoted  by  y, 
where  r  G  y  is  a  possible  state  of  y,  and  a  decision  variable  v  e  ir,  where  n  is  the  set 
of  decisions  at  each  stage.  Assume  T  decision  stages  and  let  the  variable  t  denote  the 
index  of  the  current  decision  stage.  For  a  given  stage  and  decision,  a  payoff  function 
exists  denoted  by  a(x,v,t);  and  the  state  x'  at  stage  t  +  1  is  given  by  the  function 
g( x,  v,  t ),  where 

x'  =  g(x,v,t).  (2.7) 

Let  f(x,t)  be  defined  as  the  maximum  (minimum)  value  of  /  attainable  given 
state  x  at  stage  t.  Without  loss  of  generality,  assume  the  objective  is  to  maximize 
f(xo,to),  where  x0  and  t0  are  the  initial  state  and  stage,  respectively.  Bellman’s 
equation  can  be  written  either  as  either  a  forward  or  backward  recursion.  As  a 
forward  recursion,  Bellman’s  equation  is 

f(x,t)=  max  {a:  G  y  :  a(x,v,t)  +  f[g(x,u,t),t  +  1]}  .  (2.8) 

isEn 

A  boundary  condition  is  required  to  specify  a  terminal  value  of  the  recursion.  Typ¬ 
ically,  this  is 

f(x,T)  =  0,  V  x  ex-  (2.9) 


2-6 


The  pairs  (x,t)  form  unique  nodes  in  the  dynamic  programming  state-space.  The 
key  principle  in  dynamic  programming  is  the  principle  of  optimality.  Consider  an 
optimal  path  from  node  (aq,ti)  to  node  (cc4,£4).  For  any  two  intermediate  nodes 
( x2,t2 )  and  (x3,t3)  in  this  path,  the  principle  of  optimality  states  that  the  optimal 
path  from  (x2,t2)  to  (^3,t3)  is  the  same  as  one  contained  in  the  optimal  path  from 
(aq,  h)  to  (*4,t4). 

The  primary  advantage  of  dynamic  programming  is  its  generality  as  virtually 
no  assumptions  are  imposed  on  the  underlying  functions  or  state-space.  Unfortu¬ 
nately,  the  computational  and  storage  requirements  associated  with  dynamic  pro¬ 
gramming  are  often  massive.  Dynamic  programming-based  techniques  have  been 
applied  to  both  the  one-  and  multidimensional  knapsack  problems.  We  first  consider 
one-dimensional  knapsack  applications. 

An  elementary  application  of  dynamic  programming  to  the  knapsack  problem 
is  found  in  Denardo  [5]  in  which  the  decision  v  is  to  include  one  of  the  T  items  at 
each  stage.  In  addition  a  slack  item  0  is  assumed  to  exist,  where  Co  =  0  and  a 0  =  1. 
The  state-space  is  scalar  x ,  where  x  is  the  remaining  capacity  of  the  knapsack. 
Therefore,  feasible  values  of  x  are  x  G  {0,1,2,...,  b}.  Note  that  the  stage  index  t  is 
unnecessary  in  this  particular  formulation  and  will  be  omitted.  If  item  i  is  included; 
that  is,  v  —  i,  the  payoff  function  a(x,  i )  =  q,  and  g[x,  i)  =  x  —  Wi,  i  =  1,2, ...  ,1, 
x  <  b.  The  objective  function  f(x)  is  the  maximum  value  the  knapsack  can  contain 
given  remaining  capacity  x.  A  backward  recursion  is  written  as 

f(x)  =  max  {x  >  0  :  q  +  f(x  —  wf)}  .  (2-10) 

i 

Since  no  items  can  be  inserted  into  a  knapsack  with  a  remaining  capacity  of  zero, 
the  boundary  condition  is  /( 0)  =  0,  and  the  solution  to  the  knapsack  problem  is 
f(b).  The  most  straightforward  way  to  solve  the  dynamic  program  is  by  explicit 
enumeration  in  which  each  f(x)  is  evaluated  for  all  T  +  1  possible  item  inclusions. 
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This  means  there  are  T  +  1  branchings  from  initial  state  b.  Furthermore,  each 
resulting,  nonzero  state  will  potentially  require  successive  branchings  of  size  T  +  1. 
Therefore,  the  complexity  grows  rapidly. 

Although  intuitive,  explicit  enumeration  is  extremely  inefficient.  Accelerated 
reaching  is  a  method  to  eliminate  state-space  computations.  Assume  the  items  are 
sorted  such  that  C\fa\  >  Cj/a*,  i  =  2,3 ,...,/.  Let  m(x)  be  defined  as  the  lowest 
index  item  that  can  be  included  in  an  optimal  packing  of  knapsack  size  x,  where 

m(x)  =  min  {i  :  f(i)  =  a  +  f(x  -  a*)},  (2-11) 


and  m( 0)  =  T. 

It  can  be  shown  that  there  exists  x*  =  max{x  :  m(x)  >  1};  that  is,  there  is  a 
state  x*  such  that  every  knapsack  of  size  x*  or  larger  will  include  at  least  one  item 
1  in  its  optimal  packing.  It  can  also  be  shown  that 


x*  <  (ai  —  1)H  —  x*,  (2.12) 

where  H  =  maxj  {cj},  j  =  1,2, ...,/.  Therefore,  when  in  a  state  x  >  x*  only 
the  decision  involving  the  insertion  of  item  1  needs  to  be  evaluated.  The  solution  is 
accelerated  by  continuing  to  insert  item  1  until  reaching  a  state  x  <  x*.  At  this  point, 
the  remainder  of  the  DP  must  be  evaluated  using  standard  methods.  In  practice, 
reaching  is  difficult  to  apply  for  large  and  non-integer  knapsack  problems.  So  more 
advanced  methods  have  been  developed  to  apply  dynamic  programming. 

Horowitz  and  Sahni  [11]  studied  a  number  partitioning  problem  that  is  a  special 
case  of  the  one- dimensional  KP.  They  present  a  dynamic  programming  algorithm 
that,  seeking  to  reduce  storage  requirements,  splits  the  sorted  multiset  of  weights 
into  two  multisets  having  a  difference  in  cardinality  of  0  or  1.  Each  smaller  multiset 
has  a  profit  set  associated  with  it.  A  profit-maximizing  dynamic  recursion  is  applied 
to  each  multiset,  and  the  optimal  solution  is  attained  by  searching  the  end  states  of 
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each  recursion  and  determining  an  optimal  feasible  pair.  While  this  approach  relies 
only  on  dynamic  programming,  often  dynamic  programming  is  used  in  concert  with 
other  algorithms. 

Pisinger  [20]  augments  a  class  of  algorithm  known  as  a  core  algorithm  with 
dynamic  programming.  The  core  of  a  knapsack  consists  of  those  knapsack  vari¬ 
ables  whose  optimal  values  are  different  between  the  optimal  IP  solution  and  the 
LP-relaxation.  Since  solving  the  core  is  equivalent  to  solving  the  knapsack  problem, 
the  core  is  estimated.  Once  a  core  is  estimated,  Pisinger  uses  a  dynamic,  reaching 
recursion  that  moves  bi-directionally  from  the  break  item  alternatively  inserting  and 
deleting  items.  Each  state  is  compared  to  an  upper  bound  and  fathomed  accordingly. 
Martello  et  al.  [15]  refine  Pisinger’s  algorithm  using  a  core  that  can  consist  of  non- 
consecutive  items  to  better  fill  the  knapsack.  A  similar  dynamic  program  is  applied; 
however,  the  state  space  is  controlled  actively  by  decreasing  knapsack  capacity,  solv¬ 
ing  a  surrogate  relaxation,  and  deriving  an  improved  lower  bound  through  optimal 
item  to  state-space  paring. 

The  complexity  of  the  MKP  generally  precludes  the  effectiveness  of  exact  dy¬ 
namic  programming.  Bertsimas  and  Dernir  [3]  mitigate  dynamic  programming  state 
space  requirements  by  applying  approximate  dynamic  programming  techniques  to 
the  MKP.  Instead  of  applying  the  exact  dynamic  recursion  formula,  an  approxi¬ 
mation  is  used  resulting  in  less  storage  and  computation.  Bertsimas  and  Dernir 
apply  three  different  approximations.  The  Erst  is  a  base-heuristic  that  rounds  the 
successive  LP-relaxation  solutions;  the  second  involves  a  parametric  approximation 
that  samples  the  state  space;  and  the  third  is  a  nonparametric  approximation.  The 
authors  conclude  that  the  base-heuristic  approximation  provides  the  best  solutions. 

Although  dynamic  programming  provides  a  solution  method,  it  does  not  take 
advantage  of  properties  associated  with  the  problem’s  integer  programming  formu¬ 
lation.  In  the  next  section,  we  consider  IP-based  methods. 
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2.3  IP  Formulation  Results  for  KP/MKP 


Generally,  the  most  efficient  algorithms  for  solving  knapsack  problems  involve 
the  IP  formulation.  This  is  because  information  can  be  exploited  by  relaxing  the 
integer  constraints  in  both  linear  programming  (LP),  Lagrangian,  and  surrogate 
relaxations.  An  algorithm  that  successively  uses  information  from  LP-relaxations 
is  the  branch-and-bound  algorithm  [25].  Consider  the  branch-and-bound  algorithm 
applied  to  a  knapsack  problem.  At  node  0,  the  branch-and-bound  algorithm  initially 
solves  the  LP-relaxation  of  the  KP.  The  objective  value  of  the  relaxation  forms  an 
upper  bound  for  the  optimal  IP  solution  value.  If  all  variables  have  integer  values, 
the  initial  LP-relaxation  is  optimal;  otherwise,  a  non-integer  variable  Xi  is  chosen  and 
two  new  LP-relaxations  are  created  by  imposing  an  integrality  constraint  Xi  <  \_Xi\ 
and  Xi  >  |~ay]  in  each  respective  relaxation.  This  is  branching,  and  it  forms  two 
additional  nodes. 

The  algorithm  chooses  one  of  the  relaxations  and  solves  it.  The  solution  must 
be  checked  for  feasibility  in  the  IP.  If  it  is  feasible,  and  one  or  more  non-integer 
variables  exist,  integrality  constraints  are  imposed  on  another  variable  and  branch¬ 
ing  occurs  again.  If  the  solution  is  infeasible,  further  branching  is  unnecessary  as 
all  nodes  generated  by  branching  on  an  infeasible  node  will,  themselves,  be  infeasi¬ 
ble.  This  is  called  fathoming  a  node.  Assuming  feasibility  is  maintained,  continued 
branching  eventually  yields  a  feasible,  integer  solution.  The  objective  function  value 
of  this  initial,  integer  solution  forms  a  lower  bound  for  the  objective  function  value 
of  the  optimal,  KP  solution.  Now  another  unexplored  node  is  chosen  and  branch¬ 
ing  occurs  on  it.  Because  each  node  forms  an  upper  bound  for  all  successor  nodes, 
any  node  whose  objective  function  value  is  less  than  the  value  of  the  initial  integer 
solution  need  not  be  explored  further. 

While  branch-and-bound  improves  the  lower  bound  by  successively  finding 
better  integer  solutions,  the  upper  bound  given  by  the  initial  LP-relaxation  can  be 
tightened  through  the  use  of  cutting  plane  algorithms  such  as  Gomory  cuts.  The 
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two  are  combined  in  a  branch-and-cut  algorithm  which  is  similar  to  branch-and- 
bonnd  except  that  cuts  are  generated  at  successive  LP-relaxations  and  stored  in  a 
cut  pool.  The  computations  required  to  find  an  optimal  solution  may  be  reduced 
through  application  of  preprocessing  to  eliminate  variables  and  search  heuristics  to 
determine  the  best  branching  nodes. 

Faster  branch-and-bound  algorithms  for  the  one-dimensional  KP  have  been 
developed  by  exploiting  properties  unique  to  the  knapsack  problem  to  attain  tighter 
upper-bounds,  faster  computation,  or  better  candidate  solutions.  One  of  the  first 
advances  was  the  upper-bound  developed  by  Dantzig  [14]  and  was  given  in  equation 
(2.2). 

Horowitz  and  Sahni  [11]  modified  a  previous  branch-and-bound  algorithm  mak¬ 
ing  it  a  depth-wise  search,  thereby,  reducing  computing  storage  requirements  from 
exponential  to  linear.  Based  on  the  same  profit-to-cost  sorting  of  Dantzig’s  method 
[4],  it  consists  of  forward  moves  and  backtracking.  A  forward  move  adds  the  largest 
possible  set  of  ordered,  non-included  items  into  the  current  solution;  and  backtrack¬ 
ing  removes  the  last  item  added.  Martello  and  Toth  [16]  present  a  similar  algo¬ 
rithm;  however,  the  forward  move  consists  of  two  phases  involving  building  and 
saving  the  current  solution.  This  reduces  the  number  of  backtrackings.  Of  particu¬ 
lar  significance  is  their  development  of  a  tighter  upper-bound  than  that  developed  by 
Dantzig.  For  particularly  hard  problems,  Martello  and  Toth  [14]  develop  minimum 
and  maximum  cardinality  upper  bounds  computed  by  using  a  Lagrangian-relaxation. 
A  branch-and-bound  algorithm  uses  these  upper  bounds  to  fathom  nodes. 

Core  algorithms  have  also  been  developed  that  incorporate  branch-and-bound. 
Pisinger  [19]  initially  determines  the  break  item  b ,  the  core’s  center,  through  an 
efficient  partial  sorting  algorithm.  A  greedy  algorithm  finds  an  initial  solution,  and 
a  branch-and-bound  enumerates  items  from  the  break  item  outward.  Each  time 
the  branching  grows  outside  the  core,  a  nearby  interval  of  items  undergoes  variable 
reduction  before  being  added  to  the  core. 
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In  contrast  to  one-dimensional  KPs,  the  MKP  has  very  few  specialized,  IP- 
based  procedures  [8].  Soyster  et  ah  [22]  developed  an  algorithm  suited  for  MKPs 
with  many  variables  and  few  constraints.  Based  on  partitioning  the  variables  into 
integer  and  fractional  sets,  it  successively  finds  the  optimal  integer  solution  of  a  series 
of  LP-relaxations  given  a  partial  fractional  solution.  Convergence  is  not  guaranteed, 
nor  does  it  occur  quickly  for  more  than  100  variables.  Shih  [21]  presents  a  modified 
branch-and-bound  in  which  each  variable  is  relabeled  with  m  subscripts  denoting  the 
rank  of  its  profit-to-cost  ratio  for  each  of  the  m  constraints.  The  algorithm  proceeds 
in  a  greedy  manner  by  attempting  to  load  each  constraint  to  equality.  The  minimum 
of  the  m  resulting  objective  values  is  chosen  as  the  upper  bound,  and  branching  is 
done  on  the  constraint  with  the  largest  upper  bound.  The  algorithm  was  tested  on 
thirty,  five-constraint  knapsack  problems  with  less  than  100  variables  and  performed 
reasonably  well. 

Frevillc  and  Pleateau  [6]  present  a  reduction  method,  i.e.  an  algorithm  that 
attempts  to  reduce  the  problem  size.  Surrogate  relaxation  heuristics  determine  a 
lower  bound,  and  tests  are  applied  to  eliminate  variables  and  constraints.  Frevillc 
and  Plateau  [7]  also  develop  an  effective  method  to  find  the  exact  solution  of  the 
bidimensional  knapsack  problem.  Using  the  surrogate  dual  they  show  that  an  exact 
optimal  solution  for  the  bidimensional  case  is  equivalent  to  a  search  on  the  [0, 1] 
interval.  Furthermore,  they  modify  a  previous  search  method  and  prove  optimality 
occurs  after  a  finite  number  of  iterations. 

Using  exact  algorithms  to  solve  knapsack  problems  is  often  of  greater  theoreti¬ 
cal  than  practical  interest.  Even  the  most  efficient  algorithms  can  take  an  inordinate 
amount  of  time  to  solve  larger  problems.  Heuristic  algorithms  provide  an  alternative 
and  are  discussed  in  the  next  section. 
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2-4  Potential  Heuristic  Solution  Approaches 


Knapsack-type  problems  are  ideal  for  heuristics  because  of  the  time  and  com¬ 
putation  often  involved  in  exact  solutions.  Heuristics  are  solution  methods  that  find 
optimal  or  near-optimal  solutions  but  are  not  guaranteed  to  converge  to  optimality. 
There  are  a  number  of  general  types  of  heuristic  methods  that  can  be  applied  to 
many  different  types  of  problems.  These  include  simulated  annealing,  tabu  search, 
genetic  algorithms,  etc.  Heuristics  begin  with  some  initial  point  in  the  solution  space 
and  consist  of  two  primary  phases:  performing  a  global  search  and  performing  a  lo¬ 
cal  search.  The  global  search  enables  the  heuristic  to  explore  the  solution  space  in 
order  to  avoid  settling  on  the  nearest  local  optimum.  The  local  search  allows  the 
heuristic  to  narrow  its  scope  and  seek  the  best  solution  in  a  particular  area  of  the  so¬ 
lution  space.  In  general,  a  heuristic  is  distinguished  by  the  methods  used  to  execute 
global  and  local  searches  as  well  as  the  rules  used  to  govern  transitions  between  each 
type  of  search.  A  widely  applied  heuristic  that  searches  using  probabilistic  moves  is 
simulated  annealing. 

The  simulated  annealing  heuristic  is  designed  to  mimic  the  process  of  anneal¬ 
ing  metals.  When  a  metal  is  annealed,  it  begins  in  a  molten  state.  Although  metal 
atoms  prefer  to  be  aligned,  the  high-temperature  enables  them  to  have  random  ori¬ 
entations.  Slowly  cooling  the  metal  enables  the  atoms  to  gradually  settle  into  a 
lower-energy  aligned  state;  however,  a  rapid,  cooling  schedule  can  result  in  atoms 
becoming  trapped  in  unaligned,  high-energy  states  creating  imperfections.  The  simu¬ 
lated  annealing  algorithm  behaves  similarly  in  that,  at  high  temperatures,  simulated 
annealing  allows  free  movement  throughout  the  solution  space.  It  is  effectively  a  ran¬ 
dom  search  as  suboptimal  solutions  have  a  high  probability  of  acceptance.  At  low 
temperatures,  movement  throughout  the  solution  space  is  more  restricted  and  sim¬ 
ulated  annealing  accepts  only  superior  solutions.  Therefore,  it  effectively  becomes  a 
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hill-climbing  algorithm  [18].  Consider  the  optimization  problem 

Maximize  /( x) 

Subject  to  x  G  x- 

A  simulated  annealing  heuristic  begins  at  a  high  initial  temperature,  denoted 
by  T0.  A  random  solution  x0  e  x  is  chosen  as  a  starting  point,  and  a  random  move 
is  made  to  X\  G  x,  where  X\  is  an  adjacent  or  neighboring  solution  of  x$.  Then 
the  objective  values  f(x0 )  and  f(x i)  are  compared.  Define  5  =  f(x i)  —  f(x0).  The 
neighboring  solution  can  be  accepted  under  two  scenarios:  (i)  if  S  >  0,  move  to  a?!, 

_s_ 

and  (ii)  if  5  <0,  move  to  aq  if  q  <  eTo ,  where  q  r\ j  U[0,1], 

At  this  point,  additional  solutions  can  be  explored  at  a  temperature  of  T0,  or 
the  temperature  can  be  lowered  to  Tj  according  to  some  schedule.  Let  IV  >  1  be 
defined  as  the  number  of  solutions  explored  at  each  fixed  temperature.  At  each  tem¬ 
perature  decrease,  the  heuristic  continues  from  its  last  feasible  solution  and  explores 
N  —  1  additional  solutions.  Simulated  annealing  terminates  when  the  temperature 
decreases  to  a  predetermined  terminal  temperature  Tf  <  Tq  and  N  solutions  have 
been  explored  at  Tf.  Aside  from  selection  of  the  parameters  themselves,  the  vari¬ 
ations  of  simulated  annealing  abound.  Differences  include  temperature  schedules, 
solution  selection,  termination  conditions,  and  post-processing  [18]. 

When  using  the  simulated  annealing  heuristic,  optimality  is  guaranteed  only 
under  very  specific  conditions.  Hajek  [10]  provides  a  simulated  annealing  algorithm 
that,  under  certain  conditions,  is  guaranteed  to  find  a  global  minimum  using  tem¬ 
perature  schedule  T/c  =  Cj  ln(l  +  k ),  where  C  is  the  depth  of  the  deepest  local 
minimum.  Although  of  theoretical  interest,  the  temperature  schedule  is  too  slow 
to  be  practically  useful.  In  general,  the  more  closely  a  temperature  schedule  main¬ 
tains  thermal  equilibrium  throughout,  that  is,  the  probability  distribution  of  state 
transitions  remain  close  to  their  equilibrium  distribution  at  a  given  temperature,  the 
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higher  the  quality  of  the  resulting  solution.  Triki  et  al.  [23]  provide  a  summary  of 
effective  cooling  schedules.  These  include  the  geometric  cooling  schedule  and  cooling 
schedules  that  rely  on  past  information  called  adaptive  cooling  schedules. 

The  selection  of  parameter  values  is  integral  to  the  performance  of  simulated 
annealing.  Ben-Ameur  [1]  created  an  algorithm  to  select  Tf  based  a  specified  accep¬ 
tance  probability  of  a  suboptimal  move.  The  method  requires  sampling  transitions  in 
the  state  space  and  use  of  a  recursive  formula.  More  commonly,  parameter  selection 
is  done  experimentally.  Wang  and  Wu  [24]  devised  a  six-step  method  to  experi¬ 
mentally  determine  parameter  methods  subject  to  time  constraints  using  response 
surface  methodology.  Although  theoretical  selection  of  parameters  is  not  without 
merit,  selection  of  parameters  is  often  done  experimentally  by  trying  a  variety  of 
parameter  settings  and  observing  satisfactory,  near-optimal  solutions. 

This  chapter  introduced  current  satellite  payload  selection  methodologies.  A 
general,  ten-step  process  was  discussed  to  determine  payload  specifications  from  a 
satellite’s  mission  objective.  Also  reviewed  was  a  more  specific,  quantitative  method¬ 
ology  applied  to  an  analogous  problem  of  monitoring  secondary  payloads.  The  knap¬ 
sack  problem  was  introduced  in  both  one-  and  multidimensional  cases  as  an  approach 
to  formulate  the  payload  selection  problem.  Analytical  results  for  knapsack  problems 
were  presented  as  well  as  dynamic  programming  and  integer  programming  based  so¬ 
lution  procedures.  Finally,  the  simulated  annealing  heuristic  was  introduced  as  a 
possible  solution  method  for  the  payload  selection  problem.  In  the  next  chapter,  for¬ 
mal  mathematical  models  are  presented  for  the  payload  selection  and  specification 
problem.  Four  single-satellite  models  are  presented  and  shown  to  be  equivalent  to  a 
relaxation  of  the  multi-choice,  multidimensional  knapsack  problem.  Then  one  of  the 
single-satellite  models  is  extended  to  the  multi-satellite  case.  The  chapter  concludes 
by  formulating  the  multi-satellite  model  as  both  a  dynamic  program  and  an  integer 
program. 
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3.  Mathematical  Model  Description 

In  this  chapter,  formal  mathematical  models  are  developed  for  the  satellite 
payload  selection  and  MMD  specification  problem.  The  problem  assumptions  and 
model  characteristics  are  first  described  followed  by  four  mathematical  models  de¬ 
scribing  a  single-satellite  problem.  The  final  model  is  an  extension  of  the  single¬ 
satellite  model  to  the  multiple-satellite  case.  Finally,  both  dynamic  programming 
and  integer  programming  formulations  are  provided  to  solve  the  multiple-satellite 
model. 

3.1  Model  Assumptions  and  Definitions 

Assume  a  satellite  constellation  consisting  of  S  satellites  is  observed  at  fixed 
time  epochs  with  equal  inter-inspection  time  A.  Each  satellite  bus  carries  payloads, 
and  upon  inspection  the  binary  status  (functional/not  functional)  of  each  satellite’s 
payloads  is  known  with  certainty.  It  is  also  assumed  that  at  time  n  —  0  the  mean 
mission  duration  (MMD)  specifications  of  any  payloads  on  satellites  already  in  the 
constellation  are  known  as  well  as  their  time  in  service.  At  M  predetermined  epochs 
71 1, 77.2 $:•»■■■  •  i  nMi  single-satellite  buses  will  be  loaded  with  payloads  selected  from  a 
fixed  set  of  all  payloads  and  launched  into  the  constellation.  All  selected  payloads 
must  have  a  specified,  nonzero  MMD.  It  is  assumed  that  satellite  payloads  imme¬ 
diately  enter  service  upon  launch.  Associated  with  each  satellite  payload  is  a  total 
lifetime  utility.  Utility  in  the  context  of  this  research  does  not  correspond  to  the 
strict  definition  in  utility  theory.  This  will  be  fully  discussed  later  in  the  chapter. 
The  objective  is  to  maximize  the  total  lifetime  utility  of  the  constellation. 

Each  satellite  bus  has  finite  power,  cost,  weight,  and  volume  constraints.  It  is 
assumed  that  any  set  of  specified  payloads  requiring  more  resources  than  the  satellite 
bus  can  provide  is  infeasible.  A  description  of  each  finite  resource  is  provided  in  what 
follows. 
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The  satellite’s  power  subsystem  is  responsible  for  generating  and  storing  all 
electrical  energy  required  by  the  payloads.  For  earth-orbiting  satellites,  photovoltaic 
cells  convert  solar  radiation  into  electrical  energy  which  is  then  stored  in  the  battery. 
Energy  storage  is  necessary  because  during  certain  periods  of  the  satellite’s  orbit,  the 
sun  is  eclipsed  by  the  earth  and  the  solar  cells  are  unable  to  receive  solar  radiation. 
During  eclipse,  the  battery  is  the  sole  provider  of  energy  to  the  satellite.  Therefore, 
the  battery  continuously  undergoes  charge  and  discharge  cycles  as  the  satellite  passes 
in  and  out  of  eclipse.  Because  the  duration  and  periodicity  of  eclipses  depends  on  the 
satellite’s  orbit,  the  depth  of  battery  discharge  varies  from  an  average  of  15-20%  for 
the  short,  frequent  eclipse  periods  of  low  earth  orbit  to  50%  for  the  long,  infrequent 
eclipse  periods  of  geostationary  orbit.  The  battery  cycle  life  (i.e.  the  number  of 
times  it  can  be  charged  and  discharged)  is  dependent  on  the  depth  of  discharge.  A 
greater  depth  of  discharge  reduces  the  cycle  life  of  the  battery.  Therefore,  taking  the 
expected  depth  of  discharge  into  account,  the  total  lifetime  output  of  the  battery  can 
be  estimated  by  multiplying  the  expected  energy  discharge  per  cycle  by  the  cycle 
life.  For  all  formulations,  the  power  constraint  is  equivalent  to  the  lifetime  output 
of  the  satellite’s  battery  and  is  denoted  by  P  (measured  in  Watt- Years). 

An  overall  budget  is  allotted  to  a  satellite  mission  that  usually  includes  de¬ 
velopment,  construction,  launch,  and  support  of  the  satellite.  It  is  assumed  that 
a  known  portion  of  this  total  allotment  is  assigned  specifically  to  payload  procure¬ 
ment.  This  quantity,  denoted  as  C  (measured  in  $),  will  constitute  the  total  cost 
constraint. 

The  weight  of  a  satellite  is  ultimately  constrained  by  the  capability  of  the 
selected  launch  vehicle  to  place  the  satellite  in  its  required  orbit.  The  weight  of  a 
mission-capable  satellite  at  launch  is  defined  as  the  loaded  weight  and  includes  the 
weights  of  both  propellent  and  all  satellite  subsystems.  Subtracting  the  propellent 
weight  from  the  loaded  weight  yields  the  satellite  dry  weight.  Payloads  typically 
constitute  between  15%-50%  of  satellite  dry  weight.  In  this  research,  it  is  assumed 
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that  a  launch  vehicle  is  selected  and  the  weights  of  both  propellent  and  non-payload 
subsystems  are  estimated  leaving  a  known  allowance  for  total  payload  weight  denoted 
by  W  (measured  in  lb). 

A  satellite’s  volume  is  constrained  by  the  volume  of  the  launch  vehicle’s  pay- 
load  compartment.  Additionally,  the  satellite  must  conform  to  the  compartment’s 
geometry.  Because  ensuring  geometric  feasibility  drastically  increases  the  problem 
complexity,  it  is  assumed  that  a  basic  satellite  bus  design  is  used  of  known  volume 
that  ensures  geometric  feasibility  after  payload  inclusion.  The  total  payload  volume 
constraint  V  (measured  in  ft3)  is  the  difference  between  the  volumes  of  the  launch 
vehicle’s  payload  compartment  and  the  volume  of  the  satellite  bus. 

Assume  there  exist  K  different  satellite  payload  types.  Let  L\,  L2, . . . ,  Lk 
denote  the  lifetimes  of  each  satellite  payload  type.  It  shall  be  assumed  that  these 
lifetimes  are  mutually  statistically  independent.  For  each  payload  type,  assume 
there  are  £  distinct,  nonzero  MMD  specifications.  Denote  the  MMD  of  payload  i 
on  satellite  j  by  ml,  i  =  1,2, ...  ,K,  j  =  1,2  ,...,M.  Let  0*  denote  the  set  of 
nonzero  type-i  MMD  specifications,  where  |0j|  =  £,  i  —  1,2, ...  ,K.  Additionally, 
each  payload  type  i  can  be  assigned  rrij  =  0  which  is  equivalent  to  excluding  payload 
type  i  from  satellite  j. 

A  payload  is  selected  if  it  is  assigned  a  nonzero  MMD  specification.  There¬ 
fore,  the  MMD  specification  serves  as  both  a  selection  and  specification  variable. 
The  values  rnj ,  %  =  1,  2 ,...  ,K,  j  =  1,2, ...  ,M  will  serve  as  the  decision  variables 
throughout.  For  satellite  j ,  define  raJ  =  [m\,  mJ2, . . . ,  mJK\,  j  =  1,2, ... ,  M,  as  the 
row  vector  of  MMD  specifications  for  all  payload  types  on  satellite  j.  Associated 
with  each  payload  is  a  utility  function.  Utility  in  the  classical  sense  is  a  relative  scalar 
measure  that  compares  the  probabilistic  outcomes  of  two  consequences.  Keeny  and 

Raiffa  [12]  define  utility  by  considering  a  set  of  consequences,  x\,  x2, _ ,  xn,  such 

that  their  preference  order  is  X\  -<  x2  -<...-<  xn.  If  Xj  -<  xt,  then  consequence  x, 
is  preferable  to  Xj.  A  numerical  value  scaling  ^  is  assigned  to  each  x *  such  that 
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cp i  =  0,  (fn  =  1,  and  (fi  <  (f2  <  ■  ■  ■  <  <pn,  i  —  1,  2, . . . ,  n.  Now  consider  an  action  a 
that  results  in  consequence  Xi  with  probabilities  pt .  The  utility  (p  of  a  is  calculated 
as 

<P  =  ^Pi<Pi-  (3.1) 

i 

Although  the  term  payload  utility  will  be  used  throughout  this  thesis,  it  is  not 
utility  in  the  sense  defined  by  equation  (3.1).  Instead,  payload  utility  is  analogous  to 
a  relative  value  scaling  on  M+.  The  term  utility  is  invoked  throughout  this  research 
instead  of  value  for  two  reasons:  the  computation  of  payload  utility  is  similar  to  the 
computation  of  utility  in  the  classical  sense,  and  payload  utility  is  dependent  on  a 
payload’s  relative  importance  which  is  closely  analogous  to  the  value  scaling  used  in 
the  utility  definition.  Utility  of  a  payload  type  i  is  a  function  of  three  factors: 

1.  The  relative  importance  of  payload  type  i\ 

2.  the  MMD  specification  of  payload  type  i; 

3.  the  expected  number  of  functional,  type  i  payloads  in  the  constellation. 

The  relative  importance  of  a  payload  is  only  dependent  on  its  functional  type. 
This  can  be  thought  of  as  a  measure  of  its  relative  functional  importance  to  the 
constellation’s  mission  and  is  analogous  to  a  payload’s  value.  Define  ipi  as  the  relative 
importance  of  payload  type  i,  i  —  1,2,...,  K. 

As  a  payload’s  MMD  is  increased,  it  is  expected  that  the  payload  will  operate 
longer.  However,  it  is  also  possible  that  the  shape  of  its  utility  curve  will  change  as 
its  construction  is  more  robust.  Thus,  MMD  is  included  in  the  utility  function.  The 
rationale  for  utility  dependence  on  the  number  of  functional,  like-type  payloads  is 
to  model  diminishing  marginal  returns.  For  example,  as  the  constellation  contains 
increasing  numbers  of  a  payload  type,  further  additions  of  that  type  may  add  only 
negligible  utility  to  the  constellation.  It  may  eventually  be  optimal  to  add  a  payload 
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of  another  type,  even  if  it  is  of  lower  importance.  The  random  number  of  functional 
type-i  payloads  at  time  n  is  denoted  by  a  random  variable  Qi(n)  and  its  expectation 
by  E[Qi(n)\,  i  =  1,2,.. . ,  K,  n  >  0. 

Utility  can  be  either  static  (constant  over  time)  or  dynamic  (varying  over 
time).  In  the  case  of  dynamic  utility,  it  is  assumed  that  utility  is  a  discrete-time 
stochastic  process  changing  values  only  at  inspection  times  iA,  i  =  0, 1, . . .  .  The 
static  utility  of  payload  type  i  on  satellite  j  is  denoted  by  uK^ml,  E[Qi])-,  and 
uK'ipi,  mj,  E[Qi(n)\]  n)  denotes  the  dynamic  utility  of  payload  type  %  on  satellite  j  at 
epoch  n,  i  =  1, 2, _ ,  K,  j  —  1,  2, . , ,  >  M . 

Whether  the  utility  is  static  or  dynamic,  the  total  lifetime  utility  of  a  payload 
type  %  on  satellite  j  is  represented  by  U\ (ml),  i  =  1, 2, . . . ,  K,  j  =  1,  2, ... ,  M.  For 
satellite  j,  define  the  vector 

uJ(m])  =  [Ul(m{),Ul(mJ2),...,WK^mJK)\,  j  =  1,2,  (3.2) 

The  elements  of  riJ(mJ)  are  the  total  lifetime  utilities  of  payloads  on  satellite  j 
with  MMD  specification  md,  j  =  1,2, .. . ,  M.  Each  payload  consumes  power,  cost, 
weight,  and  volume  resources.  Satellite  payloads  are,  in  general,  custom-built  items 
and  are  not  mass  produced.  The  nuances  associated  with  different  satellites  preclude 
them  having  a  standardized  design.  Estimates  are  typically  given  by  engineers  for  a 
payload’s  resource  requirements  if  designed  to  meet  a  particular  set  of  MMD  spec¬ 
ifications.  Therefore,  these  estimates  will  be  modeled  as  discrete  functions.  Define 
the  discrete  functions  Cj(m^),  Wt (ml ) ,  and  V)(ml)  that  represent  the  cost,  weight, 
and  volume,  respectively,  of  payload  type  i  with  specification  m\,  and  since  ml  =  0 
corresponds  to  excluding  payload  %  on  satellite  j,  assume  that 

C'i(O)  =  u(0)  =  Wi(o)  =  0,  i  =  l,2,...,K,  j  =  1,2, ,  M.  (3.3) 
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It  is  assumed  that  the  rate  of  a  payload’s  power  consumption  is  proportional 

to  its  utility  and  independent  of  the  payload’s  MMD  specification.  Let  Ai  be  defined 

as  the  rate  of  power  consumption  of  a  new  type-?'  payload,  i  —  1,2, ,  K.  If  utility 

is  static,  power  consumption  can  be  calculated  using  At,  i  =  1,2, ,  K.  To  allow 

for  dynamic  utility,  payload  power  must  be  calculated  by  multiplying  utility  by  a 

/  • 

scaling  factor.  Let  A-  denote  the  power  scaling  factor  of  payload  i  on  satellite  j 
where, 


A?  = 


Ai 


Ui(ipi,ml,  0;0)’ 


i  =  1,2, ...  ,K,  ]  =  1,2  ,...,M. 


(3.4) 


The  scaling  factor  is  calculated  by  multiplying  the  initial  rate  of  power  consumption 
by  l/uifyl,  mli  0;  0),  which  is  the  utility  of  payload  i  on  satellite  j  when  it  enters 
service  and  no  utility  dependence  is  present.  The  resulting  constant  maintains  the 
direct  proportionality  of  power  consumption  to  utility  when  utility  increases  or  de¬ 
creases  at  future  times;  therefore,  it  is  power  required  per  unit  of  utility. 

Total  lifetime  power  consumption  of  payload  type  i  on  satellite  j  is  defined  by 

Rj (ml),  i  —  1, 2, , . . ,  K,  j  —  1,  2, _ ,  M.  For  any  satellite  j,  define  vector  functions 

pi  (mi),  e^md),  w^(m^),  and  v^(m^)\ 


pi 

(mJ) 

■■,RK(mK)\,  J  =  1?  2,  . . 

■  ■ ,  M, 

c?\ 

(m-7) 

=  [C1(m{),C2(mJ2),. 

■■,CK(mJK)\,  j  =  1,2,.. 

■ ,  M, 

wl 

(mJ) 

■  ■  ■  ,WK{m3K)},  j  =  1,2, 

...,M, 

vl 

(md) 

=  [yi(mi),y2(mJ2),,. 

■  ■  ,VK(m3K)],  j  =  1,2,... 

,  M. 

3.2  Single- Satellite  Model 

Utility  can  also  be  characterized  as  deterministic  or  stochastic.  Any  mathe¬ 
matical  model  formulation  is  affected  by  the  underlying  characterization  of  payload 
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utility.  Interchanging  static  or  dynamic  and  deterministic  or  stochastic  leads  to  four 
possible  types  of  payload  utility: 

1.  Static  and  deterministic; 

2.  Dynamic  and  deterministic; 

3.  Static  and  stochastic; 

4.  Dynamic  and  stochastic. 

In  the  next  four  sections,  the  payload  selection  and  specification  problem  is 
formulated  for  a  single-satellite  bus  to  be  launched  into  an  empty  constellation. 
A  distinct  mathematical  programming  formulation  is  constructed  for  each  type  of 
utility.  For  notational  brevity,  the  satellite  subscript  j  will  be  suppressed,  and  the 
utility  dependence  variable  E[Qi(n )]  is  omitted  since  no  utility  dependencies  exist 
in  this  case. 

3.2.1  Static  and  Deterministic  Utility 

Consider  the  problem  of  payload  selection  and  specification  for  a  single-satellite 
bus  to  be  launched  into  an  empty  constellation.  Payload  utility  is  assumed  to  be 
static  and  deterministic;  therefore,  it  is  constant  over  time  and  known  with  certainty. 
The  total  lifetime  utility  of  payload  i  is 


Ui{mi )  =  u(^i,  mt) mi,  i  =  1,2,...,  K, 
and  the  total  lifetime  power  consumption  of  payload  i  is 


Riirrii )  =  AiU(tpi ,  mi)mu  i  =  l,2,...,K. 


(3.5) 


(3.6) 
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The  mathematical  programming  formulation  under  static  and  deterministic 


utility  is  as  follows: 


Maximize 

u(m) 

•  1 

(3.7) 

Subject  to 

p(m) 

•  1 

< 

P 

(3.8) 

c(m) 

•  1 

< 

c 

(3.9) 

v(m) 

•  1 

< 

V 

(3.10) 

w(m) 

•  1 

< 

w 

(3.11) 

mt 

G 

©iU{0},  *  =  l,2,...,/l 

(3.12) 

where  1  is  a  column  vector  of  ones. 

The  payload  selection  and  specification  problem  is  a  relaxation  of  the  MCMKP 
of  equations  (2.3)-(2.6).  The  set  of  all  items  0  =  ©i  U  @2  U  . . .  U  0 k  consists  of  each 
unique  payload  type  and  nonzero  MMD  specification  combination.  The  items  are 
partitioned  by  payload  type  into  subsets  0*,  *  =  1,2,...,  K .  The  problem  considered 
in  this  thesis  is  a  slight  relaxation  of  the  MCMKP  because  it  is  not  required  that 
one  item  of  each  type  (or  partition)  be  selected  as  in  the  MCMKP.  Payloads  can  be 
given  a  MMD  of  zero  so,  at  most  one  item  of  each  type  can  be  selected.  Power,  cost, 
weight,  and  volume  represent  multiple  dimensions  of  the  knapsack. 

Static  and  deterministic  utility  results  in  a  relatively  straightforward  formula¬ 
tion.  Because  utility  is  static  and  deterministic,  the  computation  of  total  payload 
utility  and  power  consumption  require  only  one  evaluation  of  the  utility  function 
along  with  two  and  three  multiplications,  respectively.  Therefore,  with  the  excep¬ 
tion  of  a  few  minor  calculations,  the  static  and  deterministic  case  of  payload  selection 
is  virtually  a  relaxed  MCMKP.  Next  dynamic  and  deterministic  utility  is  considered. 
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3.2.2  Dynamic  and  Deterministic  Utility 

In  this  section,  the  case  of  dynamic  and  deterministic  payload  utility  is  consid¬ 
ered.  Recall,  it  is  assumed  that  utility  is  a  discrete-time  stochastic  process.  Thus,  the 
utility  of  payload  i  remains  constant  on  all  inter-inspection  intervals.  Therefore,  the 
total  lifetime  utility  of  a  payload  is  calculated  by  summing  the  total  utility  on  each 
individual  time  interval  over  all  time  intervals  spanned  by  the  payload’s  MMD.  Let  Ki 
be  defined  as  the  number  of  time  intervals  spanned  by  payload  i ,  where  Ki  =  \rrii/  AJ , 
i  =  1,  2, . . . ,  K.  Then,  the  total  lifetime  utility  and  power  consumption  of  payload 
type  i  are  respectively  approximated  as 

Ki~  1 

Uiiymi)  =  Y,  Uj(i/;i,mi-,n)A,  (3.13) 

n= 0 

Ki- 1 

Ri{rrii)  =  Y  A'jUi^m^n) A.  (3.14) 

n=0 

Substituting  equations  (3.13)  and  (3.14)  into  vectors  u(m)  and  p{m )  leads 
to  the  same  formulation  as  equations  (3.7)-(3.12).  Therefore,  the  dynamic  and 
deterministic  utility  formulation  is  also  equivalent  to  a  relaxation  of  the  MCMKP. 
In  the  next  section,  static  and  stochastic  utility  is  considered. 

3.2.3  Static  and  Stochastic  Utility 

A  payload  with  static  and  stochastic  utility  will  operate  at  a  fixed-value  of 
utility  throughout  its  lifetime;  however,  the  exact  value  it  operates  at  is  not  known 
with  certainty.  Instead  the  possible  values  are  described  by  a  probability  distribution. 
The  expression  for  total  lifetime  utility  of  a  payload  is  identical  in  form  to  that  of 
the  static,  deterministic  case: 


Uiirrii )  =  mi)mi%  i  =  1,2,...,  K. 


(3.15) 
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However,  because  the  utility  function  is  stochastic,  it  is  advantageous  to  con¬ 


sider  the  mathematical  expectation  of  total  payload  lifetime  utility: 


E[Ui(mi)\  =  E  [ui  (ipi ,  rrii )  mi] 

=  E[ui(^i,mi)\mi ,  i  =  1,2,.. K.  (3.16) 

Similarly,  the  expected  total  lifetime  power  consumption  of  payload  i  is  given  by: 

E[Ri(mi)]  =  AiElu^i,  i  =  1,  2, ... ,  K.  (3.17) 

Now,  define  vector  functions  uE(m)  and  pE(m )  as  follows: 

uE{m)  =  [E[U1(m1)],E[U2{m2)],. . .  ,E[UK(mK)]],  (3.18) 

pE{m)  =  [E[R1(mi)},E[R2(m2)],,,.,E[RK(mK)]].  (3.19) 


The  static  and  stochastic  payload  selection  and  specification  problem  is  formulated 
by  inserting  vectors  uE  and  pE  into  equations  (3.7)  and  (3.8).  The  static  and 
stochastic  problem  formulation  is  virtually  identical  to  the  static  and  deterministic 
formulation  with  the  exception  of  the  presence  of  the  expectation  operator.  There¬ 
fore,  it  is  also  equivalent  to  the  relaxed  MCMKP.  Finally,  the  dynamic  and  stochastic 
model  is  formulated. 

3.2.4  Dynamic  and  Stochastic  Utility 

Dynamic  and  stochastic  utility  changes  in  value  over  time;  however,  its  value 
at  a  future  time  is  not  known  with  certainty  and  is  described  by  a  probability  dis¬ 
tribution.  Because  utility  is  assumed  to  be  a  discrete-time  stochastic  process,  its 
expected  value  is  constant  over  all  inter-inspection  intervals.  Therefore,  analogous 
to  the  derivation  of  the  dynamic  and  deterministic  case,  total  payload  lifetime  utility 
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and  power  consumption  for  payload  type  i  are 


EpiimM 

fvi-l 

=  E[ui{^i,mi\n)]  A,  7  =  1,2,...,  K, 

71=0 

(3.20) 

E[Ri(uii )] 

Ki-1 

=  ^  m,;  n)]  A,  i  =  1, 2, . . . ,  K, 

71= 0 

(3.21) 

respectively,  where  the  expected  utility  and  power  consumption  on  each  time  interval 
is  summed  over  all  intervals  spanned  by  the  payload  MMD. 

Substitution  of  (3.20)  and  (3.21)  into  vectors  ri^(m)  and  pE(m)  leads  to 
the  identical  formulation  of  equations  (3.7)-(3.12).  Like  the  other  three  cases,  the 
dynamic  and  stochastic  formulation  is  equivalent  to  a  relaxed  MCMKP. 

The  dynamic  and  stochastic  formulation  is  of  the  most  interest  for  the  satellite 
payload  selection  and  specification  problem  because  it  most  accurately  reflects  the 
realistic  behavior  of  satellite  payloads  in  space.  Payload  utility  changes  over  time  and 
typically  decreases  as  payload  components  fail  or  degrade  in  the  space  environment. 
Furthermore,  utility  at  any  future  time  may  not  be  known  with  certainty.  A  payload 
may  fail  immediately  and  have  zero  utility  upon  entering  service,  or  it  may  operate 
at  near-maximum  utility  throughout  its  entire  mean  mission  duration.  So,  because 
of  the  dynamic  and  stochastic  nature  of  actual  payload  utility,  the  dynamic  and 
stochastic  model  will  be  the  central  focus  of  the  remainder  of  this  thesis. 

3.3  Multi- Satellite  Extension 

Consider  an  extension  of  the  dynamic  and  stochastic  model  to  the  full-constellation 
problem  in  which  payloads  for  M  satellites  are  selected  and  specified.  Recall  that 
the  payloads  are  launched  sequentially  at  predetermined  epochs  into  a  pre-existing 
constellation  of  S  satellites.  It  is  assumed  that  both  the  launch  times  of  the  pre¬ 
existing  satellites  and  their  payload  MMDs  are  known.  Let  m\  be  defined  as  the 
remaining  MMD  of  payload  type  i  on  preexisting  satellite  j  and  be 
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defined  as  the  number  of  time  intervals  spanned  by  the  remaining  MMD  of  payload 
type  i  on  preexisting  satellite  j,  i  —  1,2 , ,K,  j  =  1, 2, . . . ,  S. 

Therefore,  U-,  the  expected  total  remaining  utility  of  payload  type  i  on  pre¬ 
existing  satellite  j  is 

Rj-l 

E[U-  (m-)]  =  j^^[wi(V’i,^,£?[Qi(n)];n)] A,  (3.22) 

rc=0 

for  *  =  1,2 j  =  1,2, ...  ,S. 

Let  uE(fhj)  be  a  vector  function  of  the  expected  total  remaining  utilities  of 
all  payloads  on  a  preexisting  satellite  j  defined  by 

=  [E[U{(m{)\,  E[Ui{m{)\, ...,  E[U3K(m3K)]],  j  —  1,2, . . . ,  S.  (3.23) 


Likewise,  for  payloads  on  the  M  satellites  to  be  launched,  the  expected  total  utility 
and  power  consumption  are 


£[bw')i 

k\-1 

=  E[ui(ipi,m3i,E[Qi(n%n) ]A, 

n= 0 

(3.24) 

EWM)} 

Kj- 1 
l 

=  Y  mhE[Qi(,n')\'in)\ A, 

n= 0 

(3.25) 

where,  n'  =  rij+n  adjusts  the  number  of  time  intervals  a  payload  has  been  in  service 
to  the  actual  time  interval  of  the  constellation. 

For  satellites  j  —  S  +  1,  S  +  2, . . . ,  S  +  M,  vectors  u3E(rn3)  and  pJE[rn3)  are 


uE(m J) 

=  [E[U{{m[)\,E[Ui(m32)\,.. 

(3.26) 

3  =  1,2,.. 

,.,M, 

pE(mJ) 

=  [^(mi)],^^)],.. 

(3.27) 

J  =  1,2,.. 

■  ■  ,M. 
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The  specification  of  M  satellites  results  in  M  power,  cost,  weight,  and  volume 
constraints  denoted  by  Pj,  Cj ,  Wj,  and  V3  for  satellites  j  —  S  +  1,  S  +  2, . . .  },S  +  M. 
The  payload  selection  and  specification  problem  of  a  satellite  constellation  with 
dynamic  and  stochastic  utility  is  as  follows: 


s 

Maximize  u3E(mr] 

3= 1 

)  •  1 

+ 

S+M 

j=S+ 1 

1 

(3.28) 

Subject  to  pJE(mr/ 

)  •  1 

< 

Pj,  j  =  1,2,.. 

• ,  M 

(3.29) 

cJ  (m?) 

)  •  1 

< 

Cj,  j  =  1,2,.. 

■ ,  M 

(3.30) 

)  •  1 

< 

Vj,  j  =  1,2,.. 

■ ,  M 

(3.31) 

w\mP) 

)  •  1 

< 

Wjt  j  =  1,2,. 

■  ■ ,  M 

(3.32) 

j 

mJi 

e 

0<  u  {0}, 

(3.33) 

*  —  1,  2,  ■  ■  ■ ,  K, 

j  =  S+l,S  +  2,. . 

.  • ,  S  +  M. 

When  the  single-satellite  case  is  extended  to  multiple-satellites,  the  effect  of 
utility  dependence  must  be  considered.  Recall  that  utility  dependence  affects  both 
the  total  lifetime  utility  and  power  requirements  of  payloads.  If  payload  utilities  are 
assumed  to  be  independent,  then  it  is  only  necessary  to  solve  the  payload  selection 
problem  once  for  each  distinct  type  of  satellite  and  apply  each  optimal  loading  to  all 
identical  satellites.  However,  if  payload  utilities  are  dependent,  the  optimal  loading 
of  each  successive  satellite  will  depend  on  previous  satellite  loadings,  ft  is  tempting 
to  conclude  that  an  exact  solution  would  consist  of  solving  a  succession  of  MCMKPs, 
one  for  each  satellite,  using  the  appropriate  resource  requirements  for  the  payloads 
available  to  each  respective  satellite.  Such  a  solution  would  maximize  the  utility 
of  each  individual  satellite  in  sequence;  however,  it  would  not  necessarily  maximize 
the  utility  of  the  overall  constellation.  In  maximizing  the  constellation’s  utility,  it 
may  be  necessary  to  assign  a  preceding  satellite  a  suboptimal  loading  such  that  a 
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subsequent  satellite  can  be  assigned  a  loading  that  results  in  an  increased  overall 
constellation  utility. 

Therefore,  due  to  utility  dependence,  the  payload  selection  problem  cannot  be 
solved  as  a  series  of  knapsack  problems.  However,  its  similarity  to  knapsack-type 
problems  suggests  two  potential  solution  methodologies.  First,  payload  assignment 
can  be  thought  of  as  a  sequential  decision  process  in  which  the  decision  to  include 
a  particular  payload  on  a  particular  satellite  is  made  at  each  stage.  Such  a  process 
lends  itself  nicely  to  a  dynamic  programming  formulation.  Second,  it  is  possible  that 
binary  variables  can  be  used  to  define  different  payload  assignments,  and  the  problem 
can  be  formulated  as  an  integer  program.  In  the  next  section  both  approaches  are 
presented. 

3.4  Dynamic  and  Integer  Programming  Formulations 

In  this  section,  the  general  constellation  payload  specification  and  selection 
problem  is  formulated  using  both  dynamic  and  integer  programming.  Both  formu¬ 
lations  provide  a  basis  for  exact  solutions  to  the  payload  selection  problem. 

3-4- 1  Dynamic  Programming  Formulation 

First  considered  is  a  dynamic  programming  approach  in  which  decisions  are 
made  in  sequential  stages.  Such  a  method  agrees  with  the  intuitive  way  one  might 
solve  this  problem.  For  example,  one  could  begin  by  specifying  payload  1  on  satellite 
1  and  then  specify  payload  2  on  satellite  1.  Once  all  satellite  1  payloads  are  speci¬ 
fied,  payloads  on  satellite  2  can  be  specified.  Continuing  in  this  manner,  one  could 
eventually  conclude  by  specifying  payload  K  on  satellite  M.  If  these  specifications 
are  made  such  that  the  total  constellation  utility  is  maximized,  the  resulting  solution 
is  equivalent  to  a  dynamic  programming  solution. 
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More  formally,  a  decision  stage  v  G  {1,2,...,  is  the  specification  of  a  unique 
satellite  and  payload  type  combination,  where  l  =  KM.  Therefore,  the  dynamic 
programming  formulation  can  be  thought  of  as  a  problem  in  which  l  payloads  are 
specified. 

Define  p1  =  [P1,  P2, . . . ,  PM],  <h  =  [Ci,  C2, . . . ,  Cm],  «i  =  [Vf,  D2, . . . ,  Vm],  and 
tui  =  [Wi,  Vh2, . . . ,  VFm]  as  vectors  of  the  initial  power,  cost,  weight  and  volume 
resources  of  satellites  1,2,...,  M.  The  stage  variable  is  v,  the  number  of  remaining, 
unspecified  payloads,  v  —  1,  2, . . . ,  i.  As  payloads  are  added,  satellite  resources  are 
consumed.  Let  pv,  cv,  vv,  and  wv  represent  the  vectors  of  remaining  power,  cost, 
weight  and  volume  resources  of  satellites  1,  2, . . . ,  M  at  stage  v,  v  —  1,  2, . . . ,  t.  The 
objective  function  f{y,pv,cv,vv,wv )  is  defined  as  the  maximum  utility  achievable 
given  v  unspecified  payloads  and  pv ,  c„ ,  vv,  and  wv  resources  remaining  on  satellites 
1,2,...,  M,  v  =  1,  2, . . . ,  l.  The  payoff  function  is  Uv(mv )  defined  as  the  expected 
total  utility  of  payload  v  with  specification  mv,  v  —  1,2,. .  .  The  expected  total 
power  consumption  of  payload  v  and  specification  mv  is  denoted  by  Rv(mv),  and  let 
Cv(mv),  Vv(mv),  and  Wv(rnv)  represent  the  cost,  volume,  and  weight,  respectively, 
of  payload  v  with  specification  rnv ,  v  —  1,2, 

Specification  of  payload  v  on  satellite  j  consumes  power,  cost,  volume,  and 
weight  resources  defined  by  vectors  pc,  cc,  wc,  and  vc,  respectively,  where 


pc  =  Rv(mv)ejr  j  =  1,  2, . . . ,  M,  v  =  1,2  (3.34) 

cc  =  Cv{mv)ej,  j  —  1,2, ... ,  M,  v  =  1, 2, . . . ,  t,  (3.35) 

wc  =  Wv(mv)ej,  j  =  1,2, ...,  M,  v  =  1,2, ...,£,  (3.36) 

vc  =  Vv(mv)ej,  j  —  1,2, ... ,  M,  v  =  1,2,  ...,£,  (3.37) 


and  ej  is  defined  as  the  jth  elementary  vector,  j  =  1,2, ... ,  M. 
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Letting  /*  represent  the  optimal  cost-to-go  function,  a  forward  dynamic  pro¬ 
gramming  recursion  is  written  as  follows: 

f*(v,pv,cv,wv,vv)  =  max{  Uv(mv )  (3.38) 

mv 

+  f*(v  -  l,p„_i  -pc,cy_i  -  Cc,wv- 1  -  -  uc)}, 

for  v  =  1,2, ...  ,T  After  the  last  stage,  there  are  no  more  payloads  to  select  and 
further  increases  in  utility  are  not  possible.  Likewise,  if  either  power,  cost,  weight, 
or  volume  resources  are  depleted,  the  addition  of  payloads  is  impossible  and  util¬ 
ity  cannot  be  increased.  Therefore,  the  dynamic  program  has  the  following  set  of 
boundary  conditions: 


/*  (O,p0>co, 

Vo,  w0)  = 

=  0, 

V 

Po, 

Co, 

Vo 

,  > 

0, 

(3.39) 

/*! 

{v,  0,  cv, 

Vv,wv )  = 

=  0, 

V 

V 

> 

o, 

V 

Cv, 

vv, 

wv 

> 

o, 

(3.40) 

/*( 

'y,pv,  o, 

VV,Wy)  = 

=  0, 

V 

V 

> 

o, 

V 

Pv , 

Vv, 

wv 

> 

0, 

(3.41) 

/*' 

[v,Pv,<k 

,  o ,wv)  = 

=  o, 

V 

V 

> 

o, 

V 

Pvr 

C-V) 

wv 

> 

0, 

(3.42) 

r 

( V,pv,C 

v,  Vv ,  0) 

=  o, 

V 

V 

> 

o, 

V 

Pv  1 

Cvi 

,  vt 

,  > 

0. 

(3.43) 

The  optimal  solution  to  the  problem  is 

f*(£,pe,ce,we,Vi).  (3.44) 

The  state-space  of  the  payload  selection  problem  grows  rapidly.  In  stage  v,  the 
payload  selection  and  specification  problem  with  M  satellites,  K  payload  types,  and 
£  nonzero  MMD  specifications  per  payload  has  the  following  state-space  size: 

States  in  stage  v  =  (£  +  1)L  (3.45) 

States  in  final  stage  =  (£  +  1)£.  (3.46) 
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It  is  clear  that  for  relatively  small  values  of  M,  K,  and  £,  the  number  of  states  can  be 
massive.  For  example,  when  M  —  2  satellites  and  K  =  5  payloads,  each  with  £  =  3 
possible  nonzero  MMD  specifications,  the  final  stage  of  the  dynamic  program  has 
(3  +  l)5'2  =  410  =  1,048,576  states  to  consider.  Figure  3.1  graphically  depicts  the 
dynamic  programming  approach  for  a  selection  and  specification  problem  in  which 
£  =  2  nonzero  MMD  specifications  are  available  for  each  payload. 


Satellite 

i 

Satellite 

i 

Satellite 

i 

i - — 

Payload  Type  Payload  Type 

1  2 

- 1  i 

Payload  Type  Payload  Type  Payload  Type 

K  1  2 

1  1 

Payload  Type  Payload  Type  Payload  Type 

K  1  2 

i 

Payload  Type 
K 

31  Final  States 

Figure  3.1  Graphical  depiction  of  a  dynamic  programming  solution  (£  =  2). 


3. 4-2  Integer  Programming  Formulation 

Now  an  integer  programming  formulation  is  considered.  Let  m']  represent  a 
payload’s  MMD  specification  on  satellite  j,  j  =  1,2,...,  M.  Define  m!  =  [m\ ,  m'2 , . . . ,  rn'M\ 
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as  the  vector  of  a  payload’s  MMD  specifications  on  satellites  1,2, ,  M .  Let, 


x 


i 

m! 


1,  if  payload  type  *  is  given  MMD  rn\  on  satellite  1, 
m'2  on  satellite  2,. .  .,m'M  on  satellite  M 


0,  otherwise 


where,  the  total  utility  associated  with  each  xlm,  is  the  sum  of  the  utilities  of  type  * 
payloads  over  all  M  satellites  at  their  respective  MMDs.  Let  ulm,  represent  the  total 
utility  derived  from  the  vector  of  specifications  m! ,  where 


M 


U, 


=  >  U^ml),  i  —  1,2, 


j= i 


K. 


(3.47) 


Furthermore,  for  xlm,,  let  the  total  consumption  of  power,  cost,  volume,  and 
weight  resources  on  satellite  j  be  denoted  by  p1^,  c^,,  v^, ,  and  w'^, ,  respectively: 


11 

Prn> 

=  ^K), 

*  =  1,2,. 

■■,K,  j 

=  1,2,. 

,.,M, 

(3.48) 

rij 

Lm' 

=  Ciim'j), 

*  =  1,2,.. 

■,K,  j 

=  1,2,.. 

• ,  M, 

(3.49) 

ij 

=  Viim'j), 

*  =  1,2,.. 

• ,  K,  j 

=  1,2,.. 

■ ,  M, 

(3.50) 

ij 

Wrrr' 

=  Wiim'j), 

,  *  =  1,2,. 

J 

'  =  1,2,. 

■  ■  ,M. 

(3.51) 
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Letting  P3 ,  C3 ,  V3 ,  and  W}  represent  the  power,  cost,  volume,  and  weight  re- 


sources  associated  with  satellite  j  the  IP  formulation  is  as  follows: 

I< 

Maximize  Y  Y 

2=1  m' 

K 

(3.52) 

Subject  to  Y  PL'^rn'  < 

i= 1 

K 

Pj,  J  =  1,2,...,M 

(3.53) 

Y  < 

2—1 

K 

Cj,  j  =  l,2 ,...,M 

(3.54) 

Yv™,x™'  - 

i= 1 

K 

Vj,  j  —  1,  2, . . . ,  M 

(3.55) 

Y  < 

2—1 

Wj,  j  —  1,2, . . . ,  M 

(3.56) 

Y x™'  = 

1,  i  —  1,2, . . . ,  K 

(3.57) 

Xm'  G 

{0,1},  i  —  1,2, . . . ,  K. 

(3.58) 

For  a  M-satellite,  K -payload  type  specihcation  problem  in  which  each  payload 
has  £  nonzero  MMDs  available,  the  number  of  binary  variables  and  constraints  is: 

Number  of  binary  variables  =  K (£  +  1)M.  (3.59) 

Number  of  constraints  =  K  +  AM.  (3.60) 

Associated  with  each  payload  type  is  the  special  ordered  set  found  in  equation  (3.57), 
and  associated  with  each  satellite  are  the  four  resource  constraint  equations  (3.53)- 
(3.56).  Although  the  number  of  binary  variables  grows  quickly  as  the  number  of 

payload  types  and  satellites  in  a  problem  increases,  it  is,  in  most  cases,  less  than  the 

number  of  states  in  the  final  stage  of  the  problem’s  corresponding  dynamic  program¬ 
ming  formulation.  Also,  unlike  the  dynamic  programming  state  growth,  the  binary 
variables  do  not  increase  substantially  when  the  number  of  payload  types  increases. 
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Superficially,  the  IP  formulation  appears  more  efficient  than  the  DP  formulation. 
However,  the  relatively  small  number  of  constraints  to  variables  is  of  some  concern. 
This  is  because  variables  increase  the  dimensionality  of  the  solution  space,  whereas, 
constraints,  while  not  affecting  the  dimensionality,  narrow  the  scope  of  the  solution 
space.  The  M  equality  constraints  represented  by  equation  (3.57)  allow  only  one 
of  K  variables  in  each  of  the  M  groups  to  be  nonzero,  thus,  greatly  reducing  the 
solution  space.  The  presence  of  these  constraints  will  offset  some  of  the  inefficiencies 
associated  with  the  formulation  having  relatively  few  constraints. 

This  chapter  began  with  a  list  of  assumptions  for  the  payload  selection  problem 
including  assumptions  for  both  the  payloads  and  satellite  buses.  The  formulation 
of  the  problem  is  dependent  on  the  four  characterizations  of  payload  utility:  static 
and  deterministic,  dynamic  and  deterministic,  static  and  stochastic,  and  dynamic 
and  stochastic.  Mathematical  models  were  developed  for  the  single-satellite  problem 
using  each  characterization  of  utility.  It  was  shown  that  each  single-satellite  model 
can  be  transformed  into  a  relaxation  of  the  multi-choice,  multidimensional  knapsack 
problem.  Since  the  dynamic  and  stochastic  model  most  closely  represents  the  realis¬ 
tic  nature  of  payload  utility,  it  was  extended  to  a  multiple-satellite  model.  To  enable 
the  multiple-satellite  model’s  solution,  both  dynamic  and  integer  programming  for¬ 
mulations  were  derived.  Solution  methods  must  now  be  applied  to  exactly  solve  these 
formulations.  In  the  next  chapter,  such  methods  are  discussed.  Additionally,  three 
prospective  heuristics  are  introduced  to  achieve  fast,  near-optimal  solutions.  Finally, 
the  performance  of  both  the  exact  and  heuristic  methods  is  evaluated  by  solving  pay- 
load  selection  problems  using  both  notional  and  randomly-generated  data. 
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4.  Numerical  Experimentation 

In  this  chapter,  the  dynamic  and  stochastic  problem  is  solved  for  notional  and 
randomly-generated  problem  instances.  In  addition  to  using  exact  methods  to  solve 
the  dynamic  and  integer  programming  formulations  of  the  payload  selection  problem, 
four  heuristics  are  devised.  The  first  two  are  based  on  an  extension  of  the  classic, 
profit-to-weight  ratio  heuristic  of  the  one-dimensional  knapsack  problem;  the  third  is 
a  greedy  heuristic  based  on  payload  utility;  and  the  fourth  is  a  simulated  annealing 
routine.  Three  small,  medium,  and  large  problems  are  solved  using  each  method  with 
a  notional  data  set.  Problem  size  is  designated  by  the  number  of  binary  variables 
in  the  resulting  IP  formulation.  Finally,  the  smallest  problem  instance  is  solved  for 
payloads  and  satellite  buses  whose  resource  requirements  and  capacities  have  been 
randomly  generated.  Results  include  both  problem  solutions  and  a  comparison  of 
the  solution  techniques’  performance. 

4-1  Exact  Solution  Methods 

A  dynamic  program  can  be  solved  using  a  variety  of  methods  many  of  which 
rely,  to  some  degree,  on  enumeration.  Pure  enumeration  is  evaluation  of  all  states  in 
the  state  space.  Although  pure  enumeration  guarantees  optimality,  it  is  very  costly  in 
terms  of  computational  requirements  and  storage.  For  cases  with  appropriate  prob¬ 
lem  structure,  approximate  dynamic  programming  may  be  used  wherein  the  cost- 
to-go  function  is  approximated  to  avoid  complete  enumeration.  The  formulation  for 
the  payload  selection  problem  assumes  no  underlying  structure  for  the  importance, 
MMD,  and  resource  requirements  of  each  payload.  This  requires  all  payloads  to  be 
examined  or,  more  specifically,  all  stages  in  the  DP  to  be  evaluated.  Therefore,  it 
is  not  possible  to  find  good  approximations  of  the  cost-to-go  function  because  no 
inferences  can  be  made  on  attributes  of  unspecified  payloads.  Additionally,  this  pre¬ 
cludes  the  development  of  bounds  to  eliminate  suboptimal  states.  For  this  reason, 
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a  pure  enumerative  routine  with  pruning  was  used  and  carried  out  via  Matlab.  At 
each  stage,  infeasible  states  are  eliminated  resulting  in  reduced  computational  and 
storage  requirements.  The  IP  formulation  in  this  research  was  solved  using  Dash 
Optimization’s  commercial  solver  Xpress.  Xpress,  like  most  industrial  solvers,  uses 
a  variety  of  methods  to  solve  an  integer  program.  Preprocessing  is  applied  to  the 
initial  problem  to  eliminate  variables,  and  heuristics  are  used  to  determine  the  best 
starting  point  for  the  branch-and-cut  algorithm,  which  uses  both  Gomory  cuts  and 
lifted  cover  inequalities  at  branch  nodes.  The  algorithm  is  strengthened  through  the 
use  of  strong  branching  in  which  branching  samples  are  ranked,  and  the  best  branch 
is  chosen. 

4-2  Greedy  and  Simulated  Annealing  Heuristics 

In  this  section,  the  heuristic  methods  are  introduced.  Two  norm-based  heuris¬ 
tics  are  presented  which  are  based  on  extending  the  one- dimensional  knapsack  prob¬ 
lem’s  profit-to-cost  ratio  heuristic  to  a  multidimensional  knapsack.  A  greedy  heuris¬ 
tic  based  on  payload  utility  is  briefly  discussed.  Also  outlined  is  an  application  of 
the  simulated  annealing  routine  to  the  payload  selection  problem. 

4-2.1  Norm-Based  Heuristics 

The  motivation  for  the  norm-based  heuristics  introduced  in  this  section  is  the 
classic  profit-to-weight  ratio  heuristic  for  the  one-dimensional  KP.  As  stated  earlier, 
the  payload  selection  problem  for  a  single  satellite  is  a  relaxation  of  the  MCMKP. 
Consider  a  single  satellite  j  and  define  the  set  of  objects  as  all  unique  pairs  ( i,mj ), 
i  =  1,  2, . . . ,  K,  j  =  1,2,...,  M;  that  is,  every  possible  payload  type  and  MMD 
combination.  Payload  weights  are  power,  cost,  weight,  and  volume  requirements.  A 
greedy  solution  to  the  multiple-satellite  payload  selection  problem  is  to  maximize  the 
utility  of  each  satellite  successively,  in  other  words,  solving  a  succession  of  MCMKPs, 
each  involving  one  of  the  M  satellite  buses.  As  discussed  earlier,  such  a  greedy 
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solution  does  not  guarantee  optimality  of  the  overall  constellation’s  utility;  however, 
all  heuristic  methods  for  the  payload  selection  problem  take  this  greedy  approach. 

For  individual  satellite  buses,  the  one-dimensional  profit-to- weight  ratio  heuris¬ 
tic  cannot  be  directly  applied  to  the  payload  selection  problem.  First,  the  satellite 
bus  is  a  four- dimensional  knapsack  instead  of  a  one-dimensional  knapsack.  One  could 
choose  either  power,  cost,  weight,  or  volume  and  apply  the  greedy  profit-to- weight 
ratio  heuristic  based  on  that  resource.  However,  failing  to  take  all  resources  into 
account  leads  to  a  heuristic  of  questionable  effectiveness.  The  second  reason  the 
classic  profit-to-weight  heuristic  cannot  be  used  is  the  additional  restriction  that  at 
most  one  payload  of  each  type  can  be  loaded  onto  a  single  satellite  bus.  Fortunately, 
this  is  easily  resolved  by  tracking  which  payload  types  have  already  been  loaded  onto 
the  satellite  bus. 

To  resolve  the  multi-dimensionality,  a  scalar  is  required  that  provides  some 
aggregate  measure  of  a  satellite’s  resource  consumption.  This  is  done  using  a  similar 
motivation  as  the  Toyota  heuristic  for  MKPs.  Dividing  each  payload’s  utility  by  an 
aggregated  scalar  provides  an  analogous  profit-to-weight  ratio.  To  motivate  such  a 
measure,  first  consider  the  p-norrn  of  a  vector  x  =  [aq,  X2,  ■  ■  ■ ,  xn] ,  where 

IMIp  =  (X1  +  X2  +  ■  •  •  +  Xn)1/P  (4-1) 

and  p  >  1.  The  p-norm  provides  a  generalized  measure  of  distance,  and  when  p  =  2, 
it  is  the  Euclidean  distance.  For  each  satellite,  there  are  A'(£+l)  payload  type/MMD 
combinations.  Define  U- ,  Rj,  Cf,  W( ,  and  V /  as  the  utility  and  power,  cost, 
weight,  and  volume  resources  used  by  payload  type/MMD  combination  i  on  satellite 
j,  respectively,  i  =  1,  2, . . . ,  K(£  +  1),  j  =  1,  2, ... ,  M.  For  payload  type/MMD 
combination  i  on  satellite  j,  define  vector  uj]  =  [R^  /  D, ,  C/  /  C} ,  V-  j  V/ ,  W-  j  IF/] , 
i  —  1,  2, . . . ,  K(£  +  1),  j  —  1,  2, . . . ,  M.  That  is,  uj]  is  the  vector  of  ratios  of  the  pay- 
load’s  resource  requirements  to  the  bus’s  resource  capacities,  i  —  1,  2, . . . ,  K(£  +  1), 
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j  =  1,2,...,  M .  Each  resource  has  different  units,  so  scaling  them  by  their  respec¬ 
tive  satellite  capacities  makes  them  comparable.  Now,  for  some  value  of  p,  define 
constant  toj  as  the  aggregate  payload  weight,  where 

4  =  Ik  lip,  *  =  1,  2,  •  ■  ■ ,  m  +  1),  j  =  1,2,...,  M.  (4.2) 

For  satellite  j,  coj  and  U-  can  be  computed  for  all  K(£  +  1)  combinations  of 
payload  type  and  MMD  specifications,  i  —  1,2, ,  K,  j  —  1,2, ,  M.  This  forms 
the  set  of  K(£  +  1)  items.  The  combinations  (or  items)  for  each  satellite  can  be 
ordered  such  that 

-  uA,  ~  ~  wjt(£+1| 

Payloads  can  then  be  inserted  into  the  satellite  greedily  from  largest  to  smallest  ra¬ 
tios.  A  formal  description  of  the  heuristic  is  as  follows: 

For  satellites  j  —  1,2, ... ,  M . 

1.  Select  value  of  p  G  M+. 

2.  Compute  tjf  and  a jf  for  each  payload/MMD  combination  i  on  satellite  j, 
i=l,2,...,K(£  +  l),j  =  l,2,...,M. 

3.  Order  combinations  such  that:  %>%>...>  Uf' (g+1) . 

4.  For  combinations  i  —  1,  2, . . .  ,  K{^  +  1):  If  combination  %  can  be  included, 
and  the  payload  type  in  combination  i  has  not  been  previously  loaded,  include 
combination  i.  Otherwise,  exclude  combination  i. 

Extensive  testing  of  the  p-norm  heuristic  for  different  values  of  p  indicated 
that  the  heuristic  was  more  accurate  for  2  <  p  <  oo.  In  particular  a  value  of 
p  —  5  is  chosen  for  all  runs  deeming  the  heuristic  the  5-norm  heuristic.  It  may  seem 
counterintuitive  that  a  non- Euclidean  norm  yields  better  results.  However,  the  effect 
of  a  larger  p- value  in  the  norm  is  that  larger  vector  components  have  more  effect 
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on  the  norm’s  value.  Since  each  component  of  u>j  represents  the  fraction  of  total 
satellite  resources  a  payload  type/MMD  combination  consumes,  the  resource  that 
consumes  the  largest  amount  of  satellite  capacity  most  strongly  affects  the  norm’s 
value.  In  other  words,  under  the  5-norm,  payloads  receive  a  larger  aggregate  measure 
of  resource  consumption  if  they  consume  a  larger  portion  of  one  resource  than  they 
would  under  the  2-norm.  This  is  the  motivation  for  the  next  heuristic. 

Although  u>j  represents  a  measure  of  resources  required  by  a  payload  type/MMD 
combination,  it  does  not  take  into  account  the  relative  scarcity  of  each  resource.  For 
example,  if  a  satellite  bus  is  severely  limited  in  its  capacity  to  produce  power,  it 
is  intuitive  that  payloads  having  a  larger  relative  power  requirement  should  have 
a  greater  effect  on  the  norm’s  value.  This  is  accomplished  by  weighting  the  norm 
elements  by  the  relative  scarcity  of  resources.  Let  ( p ,  (Jc,  Cw->  and  ( y  represent  the 
relative  scarcity  of  power,  cost,  weight,  and  volume  resources  on  satellite  j,  respec¬ 
tively,  where 

M£+d  „ . 

E  Rl 

CJP  =  - >  J  =  1,2,...,M,  (4.3) 

A? 

M£+d  „ . 

E  c? 

C 3c  =  J  —  1)  2, ... ,  M,  (4.4) 

m+ d  „  . 

E  w-f 

Civ  =  t=1w  ./  1-2 . U.  (4.5) 

K(Z+1)  „  . 

E  v? 

C3V  =  l~lA  ,  J  =  1,2,...,M.  (4.6) 

vj 
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For  payload  type/MMD  combination  i  on  satellite  j,  define  a  weighted  2-norm  a;/', 
where 


d1  =  s/cuM/Pif2  +  (uH/c^2  +  ci-(»7/»j)2  +  (WW)2  (4.7) 

for  i  =  1,  2, . . . ,  iF(£  +  1),  j  =  1,2,...,  M.  The  scalar  a;/  forms  the  basis  for  the 
next  heuristic  using  exactly  the  same  routine  as  the  first  with  the  exception  that  u/ 
is  used  in  place  of  ufl.  This  is  called  the  weighted  norm  heuristic. 

To  determine  the  effectiveness  of  the  norm-based  heuristics,  they  will  be  com¬ 
pared  to  a  strictly  greedy  heuristic  that  bases  selection  decisions  on  the  total  utility 
of  payload  type/MMD  combinations.  The  greedy  heuristic  first  computes  and  then 
sorts  the  utility  associated  with  each  payload  type/MMD  combination  in  descending 
order,  i  —  1,  2, . . . ,  K,  j  —  1,  2, . . . ,  M.  Starting  with  the  largest  utility  combination, 
it  attempts  to  include  it,  subject  to  feasibility.  It  then  attempts  to  include  each 
successive  combination  and  finishes  by  attempting  to  include  combination  K(£+  1). 
Because  it  does  not  take  payload  resource  requirements  into  account,  the  greedy 
heuristic  is  expected  to  perform  worse  than  the  norm-based  heuristics. 

4-2.2  Simulated  Annealing  Heuristic 

Simulated  annealing  is  applied  to  each  satellite  in  the  payload  selection  problem 
in  a  greedy  manner.  Like  the  two  previous  heuristics,  it  seeks  to  maximize  the  utility 
of  each  individual  satellite  successively,  so  each  satellite  is  solved  as  an  individual 
MCMKP.  Without  loss  of  generality,  the  loading  of  a  single  satellite  is  considered  in 
the  following  discussion.  Define  Xi,  i  —  1,  2, ... ,  K(£  +  1)  as: 

{1,  if  payload  type/MMD  combination  i  is  included 
0,  otherwise 
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A  solution  x  =  [x\,  x2,  ■  ■  ■ ,  xk^+i)]  is  feasible  if  x  e  y,  where  y  is  the  set  of  all 
payload  loadings  such  that: 

1.  At  most  one  of  each  payload  type  i  is  included,  i  —  1,2, ,  K. 

2.  Total  power,  cost,  weight,  and  volume  requirements  of  included  payloads  does 
not  exceed  P,  C,  W,  and  V,  respectively. 

For  any  x  E  y,  a  neighboring  solution  is  any  x'  E  x  that  differs  from  x  by  the  inclu¬ 
sion  or  exclusion  of  one  payload  type/MMD  combination,  i  =  1,2,...,  K .  Therefore, 
a  move  will  consist  of  adding  or  removing  a  payload  subject  to  feasibility.  Recall  that 
T0  is  the  initial  temperature,  and  Tf  is  the  terminal  temperature,  where  Tf  <T0.  A 
geometric  cooling  schedule  is  used  with  rate  r  e  (0, 1).  Under  a  geometric  cooling 
schedule  Tl+i  =  rlTt ,  i  =  1,2 ,...,/,  where  /  =  min {j  :  Tf  >  PTq}.  The  geometric 
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cooling  schedule  lends  itself  well  to  simulated  annealing  because  initially  tempera¬ 
ture  decreases  rapidly  from  To  avoiding  an  inordinate  amount  of  time  being  spent 
in  the  random-search  phase.  The  cooling  schedule  then  tapers  off  as  it  reaches  its 
terminal  temperature  Tf  allowing  full  exploration  of  the  local  optimum.  At  each 
temperature  level,  N  solutions  are  explored.  A  detailed  outline  of  the  heuristic  is 
presented  in  Figure  4.1. 

Selection  of  parameters  To,  Tf,  r,  and  N  is  integral  to  the  performance  of 
simulated  annealing.  Although  theoretical  selection  of  parameters  is  not  without 
merit,  an  experimental  approach  was  taken  to  determine  parameter  settings.  A 
variety  of  parameter  settings  was  tested  and  those  selected  that  gave  satisfactory, 
near-optimal  solutions.  The  simulated  annealing  parameters  used  are  found  in  Table 
4.1 


Table  4.1  Simulated  Annealing  parameter  settings. 


Parameter 

T0 

Tf 

r 

N 

Value 

500 

0.01 

0.05 

5 

4-7 


Select:  N  >  1,  r  G  (0, 1),  T0,  Tf  s.t.  0  <  T0  <  Tf 
Set:  T  =  T0,  n  —  0 
a?!  =  ranclorn(cc)  G  x 

while  T  >  Tf  do 
while  n  <  N  do 

x2  =  ranclorn(cc)  G  x 

$  =  f(x 2)  -  /(*i) 

if  6  >  0  then 

Xi  =  x2 

else  if  random(g)  U[0,1]  <  er  then 
Xi  =  x2 

end  if 
n  —  n  +  1 

end  while 
T  =  rT 
end  while 

Figure  4.1  Simulated  annealing  algorithm. 

In  simulated  annealing,  it  is  desirable  to  set  T0  large  enough  to  get  some  free 
movement  throughout  the  solution  space.  However,  a  sufficiently  high  value  may 
require  a  large  number  of  iterations  before  the  system  cools  to  Tf.  Therefore,  because 
sufficient,  free  movement  throughout  the  solution  space  cannot  always  be  practically 
ensured,  simulated  annealing  is  often  dependent  on  its  starting  solution.  A  way 
to  mitigate  this  is  to  run  multiple  replications  of  the  simulated  annealing  heuristic 
and  select  the  best  solution.  A  faster,  less  exact  simulated  annealing  heuristic  run 
over  multiple  replications  is  capable  of  providing  better  solutions  than  one  lengthy, 
computationally-expensive  routine. 
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Because  simulated  annealing,  as  it  pertains  to  the  payload  specification  prob¬ 
lem,  shows  strong  dependence  on  the  initial  solution,  replications  were  performed. 
Loadings  for  each  of  the  individual  M  satellites  were  specified  by  a  single,  simulated 
annealing  run.  However,  the  overall  M  satellite  loading  problem  was  run  multi¬ 
ple  times.  As  with  the  selection  of  parameters,  the  optimal  number  of  replications 
was  determined  through  experimentation  to  balance  optimality  with  run-time.  The 
number  of  replications  chosen  was  30. 

In  order  to  evaluate  the  performance  of  both  the  norm-based  heuristics  and 
simulated  annealing,  numerical  experiments  are  required  to  compare  the  performance 
of  the  heuristics  against  exact  solutions.  The  next  three  sections  describe  both  the 
construction  and  execution  of  these  experiments. 

4-3  Description  of  Experiment 

A  numerical  experiment  was  conducted  to  compare  the  performance  of  the 
exact  and  heuristic  solution  methods  using  notional  payload  data.  A  method  is 
deemed  exact  if,  in  theory,  it  is  guaranteed  to  yield  the  optimal  solution.  Exact 
methods  include  the  Matlab  enumeration  of  the  dynamic  programming  formulation 
and  the  Xpress  solver  solution  to  the  integer  programming  formulation.  Heuristic 
methods  include  simulated  annealing,  the  two  norm-based  heuristics,  and  the  greedy 
heuristic. 


Table  4.2  Summary  of  solution  methods. 


Exact 

Heuristic 

Pure  Enumeration 
Xpress  Solver 

Simulated  Annealing 
Weighted  2-Norm 
5-Norm 

Greedy 

Performance  is  measured  by  two  attributes:  solution  quality  and  run-time. 
Although  exact  methods  are  guaranteed  to  eventually  yield  optimality,  it  may  not 
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be  achievable  in  practice  clue  to  time  or  computational  limits.  The  design  of  the 
numerical  experiments  consists  of  two  major  parts.  The  first  is  selection  of  an  explicit 
functional  form  for  payload  utility.  The  second  is  generation  of  payload  and  satellite 
bus  data  as  well  as  the  problem  instances.  Each  is  next  discussed  in  turn. 


4-3.1  Payload  Utility  Function 

No  assumptions  are  made  regarding  the  functional  form  of  payload  utility  in 
the  derivation  of  the  models  for  payload  selection.  Both  the  exact  and  the  heuristic 
methods  will  work  with  any  arbitrary  function.  Let  indicator  function  b)  (n)  G  {0, 1} 
denote  the  functional  status  (up  or  down)  of  payload  i  on  satellite  j  at  time  n, 
i  =  1,  2, . . . ,  K ,  j  =  1,  2, . . . ,  M.  To  motivate  the  functional  form  of  utility,  condition 
payload  utility  on  b)  (n) : 

EWiWi,™3i,n)]=  Y  Elui(^i^mbn)\(l)i(n)  =  x]P{(%(n)  =  x}-  (4.8) 

ze{o,i} 


Both  E[uj  ("0j ,  'rn\ ;  n)  | (n)  =  x\  and  P{4>l(n)  =  x}  need  to  be  characterized. 
Assume  the  conditional  utility  expression  is  characterized  by 


=  0]  =  0 


(4.9) 

(4.10) 


where,  D]  ,  rn] ;  n)  is  a  deterministic  function  describing  the  utility  decay  of  pay- 
load  i  on  satellite  j  at  epoch  n,  i  =  1,  2, . . . ,  K ,  j  =  1,  2, . . . ,  M.  In  other  words,  it  is 
assumed  that  the  utility  decline  of  a  functioning  payload  is  a  deterministic  function 
and  that,  when  a  payload  fails,  it  has  zero  utility. 

The  constant  7  is  a  utility  dependence  parameter  such  that  when  7  =  0, 
there  is  no  dependence,  and  7  >  0  implies  dependence.  Both  the  enumerative  and 
heuristic  methods  require  the  calculation  of  a  payload’s  total  utility  when  it  enters 
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a  constellation,  so  total  utility  cannot  be  adjusted  for  future  payloads  added.  It  is 
possible  to  reformulate  the  IP  in  a  way  that  allows  the  total  utility  of  each  payload 
to  be  adjusted  for  the  addition  of  future  payloads;  however,  this  is  not  done  to  allow 
direct  comparison  of  the  IP  method  with  the  other  methods.  Therefore,  if  utility 
dependence  exists,  only  the  total  utilities  of  subsequent  payloads  are  diminished.  A 
value  of  7  =  0.5  was  used  throughout  the  numerical  experiments. 

In  selecting  a  functional  form  Dj (^,m-;n),  i  —  1,  2, . . . ,  K,  j  —  1,  2, . . . ,  M,  it 
is  assumed  that  the  utility  of  a  functional  payload  declines  exponentially  with  time, 
and 

Dj (V’i ,  ml ;  n)  =  r/^e~/3n/m^ ,  (3  >  0.  (4.11) 

The  constant  /3  is  a  tuning  parameter  to  adjust  the  shape  of  the  curve.  It  is  as¬ 
sumed  that  payload  survival  distributions  are  independent;  however,  no  additional 
assumptions  are  made  regarding  the  distributions,  and  it  is  not  necessary  to  have  a 
memoryless  distribution.  However,  because  satellite  payloads  are  largely  comprised 
of  electronic  components,  it  is  assumed  that  payload  survival  distributions  are  ex¬ 
ponentially  distributed.  Denote  the  survival  distribution  by 

P{(jr>(n)  =  1}  =  e-an/mS  a  >  0  (4.12) 

where,  a  is  another  tuning  parameter.  The  values  f3  =  |  In  0.5|  and  a  =  |  In  0.9|  were 
used  throughout  the  experiment.  These  values  imply  that  payloads  will  survive  to 
their  MMD  with  90%  probability,  and  upon  reaching  their  MMD,  they  will  operate 
at  50%  of  their  original  utility. 

4-3.2  Notional  Data  and  Probleml  Instances 

The  motivation  for  using  notional  data  as  opposed  to  actual  data  is  that  a 
notional  data  set  can  be  tailored  to  better  test  the  solution  methods.  For  example, 
if  a  specific  set  of  actual  data  were  used,  it  is  possible  that  power  is  the  only  limiting 
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resource.  In  this  case,  the  payload  selection  problem  would  effectively  reduce  to  a 
one-dimensional  knapsack  problem;  and  the  solution  methods  would  be  evaluated 
on  a  special,  one-dimensional  resource  case  as  opposed  to  a  more  general,  multi¬ 
dimensional  case.  The  notional  data  set  can  be  constructed  to  make  the  scarcity  of 
the  resources  competitive  thereby  avoiding  such  a  situation. 

Although  the  data  is  notional,  it  is  desirable  to  use  realistic  values  for  payload 
resource  requirements.  General  ranges  were  obtained  from  a  payload  engineer  for  the 
power,  cost,  weight,  and  volume  requirements  of  the  satellite  payloads  and  capacities 
of  the  satellite  buses.  In  practice,  a  payload’s  cost  estimate  is  directly  proportional 
to  its  weight;  therefore,  costs  and  weights  are  generally  correlated  throughout  the 
notional  data.  It  is  assumed  that  each  payload  can  be  assigned  an  MMD  of  3  years,  6 
years,  or  10  years.  The  same  set  of  MMD  choices  were  used  for  each  payload  type  to 
allow  an  easier  comparison  of  results,  but  it  is  not  necessary,  in  general,  that  each  type 
of  payload  selects  among  the  same  set  of  MMD  values.  A  payload  type  constructed 
to  a  higher  MMD  specification  will  likely  include  redundant  systems  or  more  robust 
materials  requiring  additional  cost,  weight,  and  volume  resources.  Therefore,  in  the 
notional  data,  as  MMD  is  increased,  the  cost,  weight,  and  volume  resources  consumed 
by  payloads  also  increases.  The  notional  data  used  in  this  research  are  presented  in 
Table  4.3. 

Ranges  were  also  obtained  for  the  resource  capacities  of  the  satellite  buses.  It 
is  assumed  throughout  that  all  satellite  buses  have  identical  resource  capacities,  so 
each  satellite’s  loading  can  be  compared  more  easily.  Satellite  bus  data  is  provided 
in  Table  4.4 

To  fully  explore  each  method’s  performance,  a  variety  of  problem  instances  are 
required,  both  simple  and  complex.  Problem  complexity  for  all  methods  is  dependent 
on  three  variables:  the  number  of  satellite  buses,  payload  types  available,  and  MMD 
specifications  available  to  each  payload  type.  Through  all  runs,  each  payload  type 
has  three  MMD  specifications  available  to  ensure  realism.  The  number  of  satellites 
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Tabic  4.3  Notional  satellite  payload  data. 


Type 

MMD  (Yr) 

Importance 

Power  (W) 

Cost  ($100k) 

Weight  (lb) 

Volume  (ftp) 

3 

10.0 

500 

425 

450 

15.0 

1 

6 

10.0 

500 

460 

475 

17.0 

10 

10.0 

500 

500 

500 

20.0 

3 

8.5 

475 

375 

400 

16.0 

2 

6 

8.5 

475 

405 

415 

18.0 

10 

8.5 

475 

430 

420 

19.5 

3 

7.5 

425 

410 

430 

10.0 

3 

6 

7.5 

425 

460 

480 

13.0 

10 

7.5 

425 

480 

495 

14.0 

3 

7.0 

260 

300 

230 

10.0 

4 

6 

7.0 

260 

350 

280 

13.0 

10 

7.0 

260 

370 

300 

14.0 

3 

6.0 

225 

370 

380 

13.0 

5 

6 

6.0 

225 

400 

390 

15.5 

10 

6.0 

225 

410 

395 

17.5 

3 

5.5 

300 

280 

240 

8.0 

6 

6 

5.5 

300 

320 

290 

9.0 

10 

5.5 

300 

380 

310 

12.0 

3 

5.0 

275 

150 

280 

7.0 

7 

6 

5.0 

275 

190 

350 

9.5 

10 

5.0 

275 

240 

410 

14.0 

3 

3.0 

175 

270 

225 

4.0 

8 

6 

3.0 

175 

310 

360 

5.5 

10 

3.0 

175 

335 

300 

8.0 

and  payload  types,  however,  are  varied.  Problem  size  is  measured  by  the  number  of 
binary  variables  in  the  resulting  IP  formulation.  This  measure  was  chosen  because, 
as  seen  later,  the  IP  formulation  was  the  only  exact  method  that  was  practical  for 
all  problem  sizes.  Small  problems  have  <  100  IP  variables;  medium  problems  have 
101  —  10,  000  IP  variables;  and  large  problems  have  >  10,  000  IP  variables.  Three 
instances  each  of  small,  medium,  and  large  problems  were  generated  to  provide  a 
good  sampling  over  each  range  of  problem  size. 

The  preexistence  of  satellites  in  the  constellation  at  time  n  =  0  does  not  affect 
problem  complexity.  In  all  problems,  it  is  assumed  there  are  no  preexisting  satel- 
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Tabic  4.4  Notional  satellite  bus  data. 


Power  Capacity  (W-Yr) 

Budget  ($100k) 

Weight  (lb) 

Volume  (fts) 

10500 

2500 

2500 

100 

Tabic  4.5  Payload  selection  problem  instances. 


Problem  Size 

Factor 

Instance  1 

Instance  2 

Instance  3 

No.  Satellites 

1 

2 

2 

Small 

No.  Payload  Types 

8 

4 

5 

No.  IP  Variables 

32 

64 

80 

No.  Satellites 

3 

4 

5 

Medium 

No.  Payload  Types 

4 

6 

6 

No.  IP  Variables 

256 

1536 

8192 

No.  Satellites 

6 

6 

7 

Large 

No.  Payload  Types 

6 

8 

8 

No.  IP  Variables 

24576 

32768 

131072 

lites;  therefore,  utility  dependence  is  solely  generated  by  the  payloads  selected  by 
each  method.  This  prevents  confounding  the  source  of  utility  dependency  effects. 
Satellites  are  launched  at  epochs  2,  5,  7,  9,  10,  11,  13,  and  15,  and  with  such  close 
launch  intervals,  utility  dependence  is  expected  to  occur.  All  problem  instances  be¬ 
gin  by  selecting  payloads  for  epoch  2,  the  first  launch  epoch.  Additional  satellites 
use  each  successive  launch  epoch.  Matlab  code  was  written  to  execute  the  dynamic 
programming  enumeration  and  all  heuristic  methods.  The  IP  formulation  was  solved 
using  Xpress  software.  All  computations  were  preformed  on  a  Dell  Precision  Work¬ 
station  with  a  2.66  GHz  Intel  Xeon  CPU  and  2  GB  memory.  In  the  next  section, 
numerical  results  are  presented  for  each  instance  of  small,  medium,  and  large  problem 
sizes. 


4-4  Numerical  Results  and  Summary 

This  section  contains  results  of  the  notional  payload  selection  problem  in¬ 
stances.  The  results  are  presented  from  the  smallest  problem  instance  to  the  largest. 
For  each  problem,  the  comparison  of  performance  measures  is  presented  followed  by 
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the  actual  payload  specifications  provided  by  each  method.  Simulated  annealing  is 
abbreviated  by  S.A.,  the  greedy  heuristic  is  abbreviated  by  G.H. ,  and  the  weighted 
norm  heuristic  is  abbreviated  by  W.  Norm.  First  presented  are  the  small  problems 
beginning  with  the  1  satellite,  8  payload  type  instance  in  Tables  4.6  and  4.7.  Note 
that  this  is  the  only  instance  lacking  utility  dependence  and  is  equivalent  to  solving 
a  pure  MCMKP. 


Table  4.6  Maximum  total  utility  (small  instance  1). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

228.7 

0.0 

0.5 

Enumeration 

228.7 

0.0 

751.39 

Simulated  Annealing 

219.1 

4.2 

22.25 

Greedy  Heuristic 

211.0 

7.7 

0.06 

5-Norm 

225.5 

1.4 

0.09 

Weighted  Norm 

209.6 

8.4 

0.12 

Table  4.7  Payload  specifications  (small  instance  1). 


Method 

Xpress 

Enum. 

S.A. 

G.H. 

5-Norm 

W.  Norm 

Satellite  No. 

1 

1 

1 

1 

1 

1 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

0 

0 

10 

0 

10 

Payload  3  MMD 

0 

0 

3 

0 

0 

0 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

Payload  5  MMD 

10 

10 

10 

3 

10 

0 

Payload  6  MMD 

3 

3 

0 

0 

10 

3 

Payload  7  MMD 

10 

10 

0 

0 

0 

0 

Payload  8  MMD 

0 

0 

3 

0 

3 

0 

Although  this  is  a  simple  instance,  the  enumerative  routine  took  a  compara¬ 
tively  inordinate  amount  of  time  to  find  an  exact  solution  when  compared  to  the 
IP-based,  Xpress  solver.  In  both  time  and  accuracy,  the  5-norm  performed  better 
overall  than  the  other  heuristics  coming  within  5%  of  the  optimal  solution.  Now 
consider  the  2-satcllite,  4-payload  type  result  in  Tables  4.8  and  4.9. 
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Tabic  4.8  Maximum  total  utility  (small  instance  2). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

414.1 

0.0 

0.5 

Enumeration 

414.1 

0.0 

1951.44 

Simulated  Annealing 

414.1 

0.0 

22.00 

Greedy  Heuristic 

385.6 

6.9 

0.06 

5-Norm 

382.0 

7.8 

0.09 

Weighted  Norm 

381.0 

8.0 

0.12 

Table  4.9  Payload  Specifications  (small  instance  2). 


Method 

Xpress 

Enum. 

S.A. 

G.H. 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

6 

10 

6 

10 

6 

10 

10 

10 

10 

6 

10 

6 

Payload  3  MMD 

6 

10 

6 

10 

6 

10 

0 

10 

0 

10 

0 

10 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

6 

10 

10 

10 

10 

Despite  a  doubling  of  IP  variables  from  32  to  64,  the  Xpress  solver  found  an 
optimal  solution  in  the  same  amount  of  time  as  the  previous  instance.  However, 
the  solution  time  of  the  enumerative  method  nearly  doubled.  Simulated  annealing 
also  found  the  optimal  solution  in  a  modest  amount  of  time.  The  two  norm-based 
heuristics  performed  worse  than  the  greedy  heuristic.  Presented  next  are  results  for 
the  2  satellite,  5  payload  type  problem  in  Tables  4.10  and  4.11. 


Table  4.10  Maximum  total  utility  (small  instance  3). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

437.6 

0.0 

0.7 

Enumeration* 

- 

- 

>288,0000 

Simulated  Annealing 

435.1 

0.6 

20.76 

Greedy  Heuristic 

405.2 

7.4 

0.07 

5-Norm 

430.1 

1.7 

0.14 

Weighted  Norm 

430.1 

1.7 

0.11 

For  this  instance  and  all  subsequent  instances,  complexity  of  the  enumeration 
precluded  it  from  finding  a  solution  in  a  reasonable  amount  of  time.  Meanwhile,  the 
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Table  4.11  Payload  specifications  (small  instance  3). 


Method 

Xpress 

Enum. 

S.A. 

G.H. 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

Payload  1  MMD 

10 

10 

- 

- 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

3 

10 

- 

- 

3 

10 

10 

10 

6 

10 

10 

3 

Payload  3  MMD 

3 

3 

- 

- 

6 

3 

0 

10 

0 

3 

0 

10 

Payload  4  MMD 

10 

10 

- 

- 

6 

10 

10 

0 

10 

10 

10 

10 

Payload  5  MMD 

10 

10 

- 

- 

10 

10 

3 

6 

10 

10 

3 

10 

Xpress  solver  took  only  0.2  s  longer  to  find  the  optimal  solution.  The  performance  of 
the  two  norm-based  heuristics  was  comparable,  and  they  were  closer  to  optimal  than 
the  greedy  heuristic.  Of  the  heuristics,  simulated  annealing  found  solutions  that 
were  closest  to  the  optimal  solution.  Presented  next  is  the  smallest,  medium-size 
problem  having  3  satellites  and  4  payload  types  in  Tables  4.12-4.14. 


Table  4.12  Maximum  total  utility  (medium  instance  1). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

590.0 

0.0 

0.4 

Enumeration 

- 

- 

- 

Simulated  Annealing 

568.6 

3.6 

26.77 

Greedy  Heuristic 

563.8 

4.4 

0.06 

5-Norm 

561.0 

4.9 

0.09 

Weighted  Norm 

561.0 

4.9 

0.13 

Table  4.13  Payload  specifications  (medium  instance  1). 


Method 

Xpress 

Enum. 

S.A. 

Satellite  No. 

1 

2 

3 

1 

2 

3 

1 

2 

3 

Payload  1  MMD 

10 

10 

10 

- 

- 

- 

6 

10 

10 

Payload  2  MMD 

6 

10 

10 

- 

- 

- 

6 

6 

10 

Payload  3  MMD 

6 

10 

10 

- 

- 

- 

10 

10 

10 

Payload  4  MMD 

10 

10 

10 

- 

- 

- 

10 

10 

10 

Although  the  underlying  IP  has  grown  to  256  variables,  Xpress  still  obtained  a 
solution  quickly.  Simulated  annealing  outperformed  the  other  heuristics  by  a  slight 
margin.  The  effects  of  utility  dependence  are  becoming  apparent.  In  all  solutions, 
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Tabic  4.14  Payload  specifications  (medium  instance  1). 


Method 

G.H. 

5- Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

1 

2 

3 

1 

2 

3 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

10 

10 

10 

10 

6 

10 

10 

6 

10 

Payload  3  MMD 

0 

10 

10 

0 

10 

10 

0 

10 

10 

Payload  4  MMD 

10 

6 

10 

10 

10 

10 

10 

10 

10 

satellite  3  is  assigned  all  four  payload  types  with  a  MMD  specification  of  10.  Note 
that  all  payloads  on  satellite  3  will  experience  decreased  utilities  due  to  the  number  of 
payloads  already  present  in  the  constellation.  Power  consumption  is  proportional  to 
utility,  so  each  payload  requires  less  power;  hence,  both  the  number  of  payloads  and 
their  MMD  specifications  can  be  increased  without  violating  the  power  constraint. 
Next,  consider  the  4  satellite,  6  payload  type  problem  in  Tables  4.15-4.17. 


Table  4.15  Maximum  total  utility  (medium  instance  2). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

874.5 

0.0 

85.0 

Enumeration 

- 

- 

- 

Simulated  Annealing 

850.1 

2.8 

47.90 

Greedy  Heuristic 

815.7 

6.7 

0.10 

5-Norm 

835.4 

4.5 

0.11 

Weighted  Norm 

822.7 

5.9 

0.15 

Table  4.16  Payload  specifications  (medium  instance  2). 


Method 

Xpress 

Enurn. 

S.A. 

Satellite  No. 

1 

2 

3 

4 

1 

2 

3 

4 

1 

2 

3 

4 

Payload  1  MMD 

6 

10 

10 

10 

- 

- 

- 

- 

10 

3 

3 

10 

Payload  2  MMD 

6 

3 

10 

10 

- 

- 

- 

- 

3 

10 

10 

10 

Payload  3  MMD 

3 

6 

0 

10 

- 

- 

- 

- 

0 

6 

10 

3 

Payload  4  MMD 

10 

10 

10 

10 

- 

- 

- 

- 

10 

6 

10 

10 

Payload  5  MMD 

10 

10 

10 

10 

- 

- 

- 

- 

6 

6 

10 

10 

Payload  6  MMD 

3 

6 

10 

0 

- 

- 

- 

- 

10 

10 

6 

10 
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Tabic  4.17  Payload  specifications  (medium  instance  2). 


Method 

G.H. 

5- Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

4 

1 

2 

3 

4 

1 

2 

3 

4 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

10 

10 

10 

10 

0 

10 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

0 

10 

10 

10 

0 

0 

10 

10 

0 

10 

0 

10 

Payload  4  MMD 

10 

0 

10 

10 

10 

10 

10 

10 

10 

0 

10 

10 

Payload  5  MMD 

3 

6 

6 

0 

10 

10 

0 

0 

0 

0 

10 

0 

Payload  6  MMD 

0 

0 

0 

10 

10 

3 

3 

10 

3 

6 

10 

10 

Finding  an  exact  solution  requires  substantially  more  time  for  the  IP  solver 
than  in  the  previous  instance.  Simulated  annealing  found  a  near-optimal  solution  in 
roughly  half  the  time.  Of  the  remaining  heuristics,  only  the  5-norm  found  a  solution 
within  5%  of  optimal.  The  largest,  medium-sized  problem  is  presented  next  in  Tables 
4.18-4.21.  This  is  the  first  instance  in  which  Xpress  fails  to  find  a  provable,  optimal 
solution.  Therefore,  all  methods  are  compared  against  the  best,  integer  solution 
obtained  as  well  as  the  lowest,  upper  bound  on  the  objective  function  value.  This 
bound  is  determined  through  cuts  generated  by  Xpress. 


Table  4.18  Maximum  total  utility  (medium  instance  3). 


Method 

Total  Utility 

%  Diff.  Best 
Solution 

%  Diff.  Upper 
Bound 

Time  (s) 

Xpress  Solver 

1088.8 

0.0 

0.2 

1404.5 

Enumeration 

- 

- 

- 

- 

Simulated  Annealing 

1039.4 

4.5 

4.7 

109.13 

Greedy  Heuristic 

1033.0 

5.1 

5.3 

0.11 

5-Norm 

1054.3 

3.2 

3.4 

0.13 

Weighted  Norm 

1037.1 

4.7 

5.0 

0.17 

The  simulated  annealing  solution  is  within  5%  of  the  best  integer  solution 
and  upper  bound;  however,  the  5-norm  came  the  closest  to  optimality  of  all  the 
heuristics.  The  largest  problems  are  next  presented  beginning  with  the  6  satellite,  6 
payload  type  instance  found  in  Tables  4.22-4.25.  Xpress  managed  to  find  a  solution, 
but  it  required  a  substantial  amount  of  time.  Although  simulated  annealing  found 
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Tabic  4.19  Payload  specifications  (medium  instance  3). 


Method 

Xpress 

Enum. 

Satellite  No. 

1 

2 

3 

4 

5 

1 

2 

3 

4 

5 

Payload  1  MMD 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

Payload  2  MMD 

3 

0 

0 

0 

10 

- 

- 

- 

- 

- 

Payload  3  MMD 

0 

3 

6 

10 

0 

- 

- 

- 

- 

- 

Payload  4  MMD 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

Payload  5  MMD 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

Payload  6  MMD 

3 

10 

10 

10 

10 

- 

- 

- 

- 

- 

Payload  7  MMD 

3 

6 

10 

10 

10 

- 

- 

- 

- 

- 

Payload  8  MMD 

0 

0 

0 

0 

0 

- 

- 

- 

- 

- 

Tabic  4.20  Payload  specifications  (medium  instance  3). 


Method 

S.A. 

G.H. 

Satellite  No. 

1 

2 

3 

4 

5 

1 

2 

3 

4 

5 

Payload  1  MMD 

6 

0 

10 

10 

6 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

10 

0 

6 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

6 

6 

3 

10 

10 

0 

10 

10 

0 

10 

Payload  4  MMD 

10 

10 

6 

10 

6 

10 

0 

10 

10 

10 

Payload  5  MMD 

6 

6 

10 

10 

0 

3 

6 

0 

10 

10 

Payload  6  MMD 

6 

3 

10 

0 

10 

0 

0 

0 

10 

0 

Payload  7  MMD 

0 

6 

10 

3 

10 

0 

0 

6 

3 

6 

Payload  8  MMD 

3 

3 

0 

0 

0 

0 

0 

0 

0 

0 

Table  4.21  Payload  specifications  (medium  instance  3). 


Method 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

4 

5 

1 

2 

3 

4 

5 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

10 

0 

10 

10 

10 

0 

10 

10 

10 

Payload  3  MMD 

0 

0 

10 

10 

0 

0 

10 

0 

10 

0 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

0 

10 

10 

Payload  5  MMD 

10 

0 

10 

10 

0 

0 

0 

10 

10 

0 

Payload  6  MMD 

10 

0 

10 

0 

10 

3 

10 

6 

0 

10 

Payload  7  MMD 

0 

10 

0 

0 

10 

0 

0 

10 

0 

10 

Payload  8  MMD 

3 

0 

0 

0 

10 

0 

3 

0 

0 

10 
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a  solution  closest  to  optimal,  the  difference  in  solution  quality  of  the  heuristics  is 
small.  Consider  the  6  satellite,  8  payload  problem  in  Tables  4.26-4.29.  In  this 
instance,  simulated  annealing  performed  worse  than  all  other  heuristics  including 
the  greedy  heuristic.  The  5-norm  performed  best.  In  the  final  instance  involving  7 
satellites  and  8  payloads  in  Tables  4.30-4.33,  the  5-norm,  once  again,  had  the  most 
optimal  solution  followed  by  the  weighted  norm. 


Table  4.22  Maximum  total  utility  (large  instance  1). 


Method 

Total  Utility 

%  Diff.  Optimal 

Time  (s) 

Xpress  Solver 

1248.0 

0.0 

5435.7 

Enumeration 

- 

- 

- 

Simulated  Annealing 

1191.4 

4.5 

71.02 

Greedy  Heuristic 

1186.9 

4.9 

0.11 

5-Norm 

1188.4 

4.8 

0.13 

Weighted  Norm 

1180.4 

5.4 

0.25 

Table  4.23  Payload  specifications  (large  instance  1). 


Method 

Xpress 

Enum. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

6 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  2  MMD 

6 

6 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  3  MMD 

3 

3 

3 

3 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  4  MMD 

10 

10 

10 

10 

3 

10 

- 

- 

- 

- 

- 

- 

Payload  5  MMD 

6 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  6  MMD 

6 

6 

6 

10 

10 

3 

- 

- 

- 

- 

- 

- 

Several  conclusions  can  be  drawn  from  the  results  of  these  notional  problem 
instances.  Of  the  exact  solution  methods,  the  IP-based  Xpress  solver  clearly  worked 
best.  Although  it  was  unable  solve  three  of  the  larger  problems  to  optimality,  it 
placed  relatively  tight  bounds  on  the  optimal  objective  function  value  and  provided 
a  near-optimal  integer  solution.  The  enumerative  method  took  inordinately  long, 
even  for  small  problems  and  if  terminated  early,  provided  no  solution.  Of  the  heuris¬ 
tic  methods,  the  simulated  annealing  routine  and  the  5-norm  provided  solutions 
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Table  4.24  Payload  specifications  (large  instance  1). 


Method 

S.A. 

G.H. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

6 

10 

6 

10 

10 

6 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

3 

10 

10 

6 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

3 

6 

10 

10 

10 

10 

0 

10 

10 

10 

10 

10 

Payload  4  MMD 

6 

6 

6 

10 

6 

10 

10 

0 

10 

10 

10 

3 

Payload  5  MMD 

10 

3 

3 

10 

10 

10 

3 

6 

6 

0 

10 

10 

Payload  6  MMD 

10 

0 

10 

6 

3 

6 

0 

0 

0 

10 

3 

10 

Tabic  4.25  Payload  specifications  (large  instance  1). 


Method 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

0 

0 

10 

10 

3 

10 

0 

10 

0 

10 

10 

10 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

0 

10 

10 

10 

10 

Payload  5  MMD 

10 

10 

0 

0 

10 

0 

0 

0 

10 

0 

10 

0 

Payload  6  MMD 

10 

3 

3 

10 

10 

10 

3 

6 

10 

10 

3 

10 

Table  4.26  Maximum  total  utility  (large  instance  2). 


Method 

Total  Utility 

%  Diff.  Best 
Solution 

%  Diff.  Upper 
Bound 

Time  (s) 

Xpress  Solver 

1286.3 

0.0 

1.0 

3472.1 

Enumeration 

- 

- 

- 

- 

Simulated  Annealing 

1219.6 

5.2 

6.1 

131.81 

Greedy  Heuristic 

1232.3 

4.2 

5.1 

0.12 

5-Norm 

1250.3 

2.8 

3.7 

0.16 

Weighted  Norm 

1234.6 

4.0 

4.9 

0.18 

that  were  consistently  close  to  optimality.  Simulated  annealing  is  the  most  time 
consuming,  yet  the  most  reliable  heuristic.  On  all  but  one  instance,  in  which  the 
sub-optimality  was  5.2%,  simulated  annealing  came  within  5%  of  the  optimal  or  best 
integer  solution. 
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Table  4.27  Payload  specifications  (large  instance  2). 


Method 

Xpress 

Enum. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

6 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  2  MMD 

0 

3 

3 

0 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  3  MMD 

6 

0 

10 

10 

10 

0 

- 

- 

- 

- 

- 

- 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  5  MMD 

10 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

Payload  6  MMD 

6 

6 

6 

10 

3 

10 

- 

- 

- 

- 

- 

- 

Payload  7  MMD 

3 

6 

0 

0 

0 

10 

- 

- 

- 

- 

- 

- 

Payload  8  MMD 

0 

3 

0 

10 

0 

0 

- 

- 

- 

- 

- 

- 

Table  4.28  Payload  specifications  (large  instance  2). 


Method 

S.A. 

G.H. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

3 

6 

10 

10 

6 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

6 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

10 

0 

6 

0 

10 

10 

0 

10 

10 

0 

10 

10 

Payload  4  MMD 

3 

6 

0 

10 

10 

3 

10 

0 

10 

10 

10 

10 

Payload  5  MMD 

10 

10 

10 

10 

0 

0 

3 

6 

0 

10 

10 

0 

Payload  6  MMD 

10 

3 

0 

3 

0 

10 

0 

0 

0 

10 

0 

10 

Payload  7  MMD 

3 

10 

3 

3 

10 

10 

0 

0 

6 

3 

6 

10 

Payload  8  MMD 

0 

3 

6 

10 

6 

0 

0 

0 

0 

0 

0 

0 

Table  4.29  Payload  specifications  (large  instance  2). 


Method 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

10 

0 

10 

10 

10 

10 

0 

10 

10 

10 

10 

Payload  3  MMD 

0 

0 

10 

10 

0 

10 

0 

10 

0 

10 

0 

10 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

0 

10 

10 

10 

Payload  5  MMD 

10 

0 

10 

10 

0 

0 

0 

0 

10 

10 

0 

0 

Payload  6  MMD 

10 

0 

10 

0 

10 

10 

3 

10 

6 

0 

10 

10 

Payload  7  MMD 

0 

10 

0 

0 

10 

10 

0 

0 

10 

0 

10 

10 

Payload  8  MMD 

3 

0 

0 

0 

10 

0 

0 

3 

0 

0 

10 

0 
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Table  4.30  Maximum  total  utility  (large  instance  3). 


Method 

Total  Utility 

%  Diff.  Best 
Solution 

%  Diff.  Upper 
Bound 

Time  (s) 

Xpress  Solver 

1469.7 

0.0 

2.8 

7012.7 

Enumeration 

- 

- 

- 

- 

Simulated  Annealing 

1415.6 

3.7 

6.3 

151.88 

Greedy  Heuristic 

1420.4 

3.4 

6.0 

0.14 

5-Norm 

1441.3 

1.9 

4.6 

0.16 

Weighted  Norm 

1426.5 

2.9 

5.6 

0.21 

Tabic  4.31  Payload  specifications  (large  instance  3). 


Method 

Xpress 

Enum. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

7 

1 

2 

3 

4 

5 

6 

7 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  2  MMD 

3 

3 

3 

0 

10 

10 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  3  MMD 

3 

0 

6 

10 

0 

10 

0 

- 

- 

- 

- 

- 

- 

- 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  5  MMD 

6 

10 

10 

10 

10 

10 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  6  MMD 

3 

3 

10 

10 

10 

0 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  7  MMD 

0 

6 

0 

0 

0 

6 

10 

- 

- 

- 

- 

- 

- 

- 

Payload  8  MMD 

0 

6 

0 

10 

3 

0 

0 

- 

- 

- 

- 

- 

- 

- 

Table  4.32  Payload  specifications  (large  instance  3). 


Method 

S.A. 

G.H. 

Satellite  No. 

1 

2 

3 

4 

5 

6 

7 

1 

2 

3 

4 

5 

6 

7 

Payload  1  MMD 

10 

6 

0 

10 

6 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

6 

10 

0 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  3  MMD 

0 

10 

10 

6 

10 

0 

0 

0 

10 

10 

0 

10 

10 

10 

Payload  4  MMD 

10 

3 

3 

6 

0 

10 

10 

10 

0 

10 

10 

10 

10 

10 

Payload  5  MMD 

6 

0 

6 

10 

10 

10 

10 

3 

6 

0 

10 

10 

0 

10 

Payload  6  MMD 

6 

10 

10 

10 

10 

0 

10 

0 

0 

0 

10 

0 

10 

0 

Payload  7  MMD 

10 

3 

3 

10 

0 

10 

10 

0 

0 

6 

3 

6 

10 

6 

Payload  8  MMD 

0 

0 

10 

0 

10 

6 

0 

0 

0 

0 

0 

0 

0 

0 

For  the  heuristic  methods  it  is  important  to  determine  whether  the  observa¬ 
tions  made  regarding  each  heuristic  method  hold  true  in  general.  For  any  heuristic, 
it  is  likely  that  it  will  perform  well  for  certain  payload  data  sets  and  poorly  for  oth- 
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Table  4.33  Payload  specifications  (large  instance  3). 


Method 

5-Norm 

W.  Norm 

Satellite  No. 

1 

2 

3 

4 

5 

6 

7 

1 

2 

3 

4 

5 

6 

7 

Payload  1  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

10 

Payload  2  MMD 

0 

10 

0 

10 

10 

10 

10 

10 

0 

10 

10 

10 

10 

10 

Payload  3  MMD 

0 

0 

10 

10 

0 

10 

10 

0 

10 

0 

10 

0 

10 

10 

Payload  4  MMD 

10 

10 

10 

10 

10 

10 

10 

10 

10 

0 

10 

10 

10 

10 

Payload  5  MMD 

10 

0 

10 

10 

0 

0 

10 

0 

0 

10 

10 

0 

0 

10 

Payload  6  MMD 

10 

0 

10 

0 

10 

10 

0 

3 

10 

6 

0 

10 

10 

0 

Payload  7  MMD 

0 

10 

0 

0 

10 

10 

6 

0 

0 

10 

0 

10 

10 

6 

Payload  8  MMD 

3 

0 

0 

0 

10 

0 

0 

0 

3 

0 

0 

10 

0 

0 

ers.  Therefore,  it  is  desirable  to  test  each  heuristic  on  a  large  number  of  randomly 
generated  data  sets  in  order  to  draw  stronger  conclusions  regarding  the  performance 
of  each  method.  This  is  done  in  the  next  section. 


4-5  Random  Problem  Instances 

In  this  section,  randomly  generated  problem  instances  are  solved  by  each 
heuristic,  and  their  solutions  are  compared  to  optimal  solutions.  Because  the  over¬ 
head  associated  with  conditioning  data  for  Xpress  is  time-consuming,  the  enumer¬ 
ation  routine  was  used  to  exactly  solve  each  problem  instance.  In  order  to  keep 
the  time  required  for  an  exact  solution  manageable,  only  the  1  satellite,  8  payload 
problem  was  solved.  In  addition  to  payloads  having  randomly-generated  resource 
requirements,  the  satellite  bus  has  randomly-generated  resource  capacities.  Two 
separate  methods  were  used  to  generate  data,  and  100  instances  of  each  method’s 
problems  were  solved. 

The  motivation  of  the  first  method  is  to  generate  a  variety  of  random  problems 
within  ranges  of  expected  payload  resource  requirements  and  satellite  bus  capacities. 
The  method,  denoted  Ml ,  uses  the  same  values  of  payload  importance  as  the  notional 
data  in  Table  4.3.  All  payload  resource  requirements  and  satellite  bus  capacities  were 
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generated  using  a  uniform  distribution.  Payloads  of  the  same  type  are  assigned  the 
same  power  requirement;  however,  unlike  the  notional  data,  there  is  no  correlation 
between  payload  MMD  and  randomly-generated  resource  consumption.  Tables  4.34 
and  4.35  show  the  ranges  over  which  the  payload  and  satellite  bus  data  was  generated: 


Table  4.34  Ranges  of  payload  requirements  for  Ml  problem  instances. 


Resource 

Power  (W) 

Cost  ($100k) 

Weight  (lb) 

Volume  (fta) 

Range 

100-500 

100-500 

100-500 

3-20 

Table  4.35  Ranges  of  bus  capacities  for  Ml  problem  instances. 


Resource 

Power  (W-Yr) 

Cost  ($100k) 

Weight  (lb) 

Volume  (fta) 

Range 

7500-12500 

2000-3000 

2000-3000 

75-125 

The  second  method,  denoted  M2,  is  designed  to  test  the  heuristics  over  a 
wider-range  of  problem  instances  than  Ml.  Instead  of  using  the  notional  vales, 
payload  importance  values  are  randomly-generated  between  0  and  10.  Data  in  M2  is 
generated  using  a  uniform  distribution  on  ranges  specified  in  Tables  4.36  and  4.37. 


Table  4.36  Ranges  of  payload  requirements  for  M2  problem  instances. 


Resource 

Power  (W) 

Cost  ($100k) 

Weight  (lb) 

Volume  (fta) 

Range 

0-500 

0-500 

0-500 

0-20 

Table  4.37  Ranges  of  bus  capacities  for  M2  problem  instances. 


Resource 

Power  (W-Yr) 

Cost  ($100k) 

Weight  (lb) 

Volume  (fta) 

Range 

0-10500 

0-2500 

0-2500 

0-100 

To  execute  the  experiment,  100  problem  instances  of  the  1  satellite,  8  payload 
type  problem  were  exactly  solved  via  DP  enumerations  using  Ml  and  M2  randomly- 
generated  data.  Then  each  heuristic  was  applied  to  the  problems,  and  the  resulting 
solutions  were  compared  to  the  exact  solutions.  The  heuristic  performance  on  the 
Ml  generated  problem  instances  are  in  Tables  4.38  and  4.39,  and  the  results  of  the 
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M2  generated  problem  instances  are  in  Tables  4.40  and  4.41.  Table  4.42  shows  how 
many  Ml  and  M2  problem  solutions  were  beyond  5%  of  optimality. 

The  mean  differences  from  optimality  of  the  heuristics’  solutions  were  relatively 
close  in  the  Ml  problem  set.  Simulated  annealing  has  the  lowest  mean  followed  by 
the  5-norm  with  a  slightly  higher  mean.  The  lower  standard  deviation  of  simulated 
annealing  indicates  that  it  provides  near-optimal  solutions  more  consistently  than 
do  the  other  heuristics.  In  the  Ml  problem  instances,  the  weighted  norm  performed 
roughly  equivalently  to  the  greedy  heuristic,  and  no  clear  advantage  was  attained 
in  using  the  weighted  norm.  The  performance  results  of  the  M2  problem  instances 
provide  a  starker  contrast  of  solution  quality  between  the  different  heuristics.  Sim¬ 
ulated  annealing  has  a  much  lower  mean  difference  from  optimality  than  any  other 
heuristic.  The  5-norm  once  again  performs  better  than  the  weighted  norm;  how¬ 
ever,  unlike  the  Ml  problem  set,  the  weighted  norm  performs  better  than  the  greedy 
heuristic. 


Table  4.38  Heuristic  solution  quality  (100  replications  of  Ml  random  data). 


Method 

%  Mean  Diff. 

%  Std.  Dev. 

%  Max  Diff. 

%  Median  Diff. 

Simulated  Annealing 

2.4931 

2.0204 

7.7330 

2.1342 

Greedy  Heuristic 

5.0545 

5.1807 

24.9068 

4.3784 

5- Norm 

3.0781 

3.5238 

15.6261 

1.8375 

Weighted  Norm 

4.9201 

5.4404 

26.9841 

3.3097 

Table  4.39  Heuristic  solution  time  (100  replications  Ml  random  data). 


Method 

Simulated  Annealing 

Greedy  Heuristic 

5-Norm 

Weighted  Norm 

Mean  Time  (s) 

20.87 

0.02 

0.02 

0.02 

A  number  of  conclusions  can  be  drawn  from  the  results  of  the  notional  and 
random  problem  instances.  Clearly,  the  exact  method  with  the  best  performance  is 
applying  the  Xpress  solver  to  the  IP  formulation.  The  DP  enumeration  routine  could 
only  solve  the  two  smallest  notional  problems  in  a  reasonable  amount  of  time.  Xpress 
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Tabic  4.40  Heuristic  solution  quality  (100  replications  of  M2  random  data). 


Method 

%  Mean  Diff. 

%  Std.  Dev. 

%  Max  Diff. 

%  Median  Diff. 

Simulated  Annealing 

0.2013 

0.8337 

5.3157 

0.0000 

Greedy  Heuristic 

10.2515 

14.1882 

74.8440 

4.0379 

5- Norm 

3.4383 

5.7982 

26.8590 

0.0000 

Weighted  Norm 

4.4508 

8.5532 

46.8890 

0.0000 

Table  4.41  Heuristic  solution  time  (100  replications  M2  random  data). 


Method 

Simulated  Annealing 

Greedy  Heuristic 

5-Norm 

Weighted  Norm 

Mean  Time  (s) 

24.26 

0.02 

0.02 

0.03 

Table  4.42  Number  of  solutions  beyond  5%  optimality  (100  replications  Ml  and  M2). 


Method 

Simulated  Annealing 

Greedy  Heuristic 

5-Norm 

Weighted  Norm 

Ml 

11 

44 

24 

40 

M2 

1 

49 

25 

25 

solved  all  but  three  of  the  notional  problem  instances,  and  for  the  three  it  failed  to 
solve,  the  best  integer  solutions  were  within  3%  of  optimality.  The  best  heuristic 
based  on  solution  quality  and  consistency  is  simulated  annealing.  The  closeness  to 
optimality  of  the  5-norm  solutions  often  rivals  those  of  simulated  annealing  but,  the 
5-norm  fails  to  attain  near-optimal  solutions  consistently.  However,  the  5-norm  ran 
on  the  order  of  1,000-10,000  times  faster  than  simulated  annealing,  and  it  can  be 
argued  that,  taking  time  into  consideration,  the  5-norm  heuristic  has  the  best  overall 
performance.  However,  practically  considering  the  lengthy  time  frame  associated 
with  satellite  construction  and  the  monetary  value  of  the  resources  at  stake,  the  time 
required  by  simulated  annealing  (~3  minutes  for  the  largest  problem)  is  negligible 
compared  to  the  value  gained  attaining  a  nearly-optimal  solution. 

Several  factors  must  be  considered  when  deciding  between  the  use  of  the  Xpress 
IP  solver  or  simulated  annealing.  Most  importantly  is  what  level  of  sub-optimality 
can  be  tolerated.  Although  this  research  provides  the  means  by  which  to  select  and 
specify  payloads,  the  solution  is  just  a  starting  point  in  the  process  of  actually  con- 
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structing  a  satellite.  Modifications  to  the  solution  will  likely  be  required  because  of 
geometric  or  thermal  considerations.  Additionally,  the  mission  objectives  or  prior¬ 
ities  of  the  satellite  constellation  itself  may  change.  Therefore,  a  5%  difference  in 
optimality  between  the  exact  and  simulated  annealing  solutions  may  not  be  very 
significant  with  respect  to  the  other  design  factors.  Additionally,  because  of  the  size 
of  the  IP  formulation,  a  costly  commercial  solver  like  Xpress  or  CPLEX  is  almost 
certainly  required.  In  contrast,  the  simulated  annealing  routine  is  well- documented 
and  coding  it  in  virtually  any  computer  language  requires  very  little  specialized  or 
proprietary  knowledge.  In  any  case,  the  tradeoff  between  a  solution’s  quality  and  its 
complexity  must  be  assessed  in  determining  which  technique  to  apply. 
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5.  Conclusions  and  Future  Research 


Nations  spend  an  enormous  amount  of  financial  resources  acquiring,  launch¬ 
ing,  and  operating  satellites.  Selecting  and  specifying  satellite  payloads  is  a  process 
that  can  benefit  tremendously  from  the  development  of  methodologies  for  more  effec¬ 
tive  resource  allocation.  Payloads  are  the  mission  critical  components  of  a  satellite 
and  require  power,  cost,  weight,  and  volume  resources,  of  which,  the  satellite  bus 
can  only  provide  a  limited  amount.  A  critical  specification  of  a  satellite  payload  is 
its  mean  mission  duration  (MMD),  which  denotes  the  payload’s  lifetime  for  mission 
planning  purposes.  Increasing  a  payload’s  MMD  specification  requires  additional 
components  and  materials  increasing  the  payload’s  cost,  weight,  and  volume.  Cur¬ 
rent  payload  selection  and  specification  methodologies  are  either  general  and  qual¬ 
itative  or  only  applicable  to  a  specific  type  of  satellite  constellation.  This  research 
has  developed  a  general,  quantitative  methodology  to  select  and  specify  satellite 
payloads  that  can  be  applied  to  virtually  any  satellite  constellation. 

The  payload  selection  and  specification  problem  was  shown  to  be  similar  to 
a  class  of  mathematical  programming  problems  known  as  knapsack  problems.  To 
motivate  a  methodology  for  the  payload  selection  problem,  a  thorough  review  of  the 
literature  on  knapsack  problems  was  presented  to  include  well-known  results  and 
a  variety  of  solution  techniques.  Solution  techniques  are  generally  based  on  two 
formulations  of  the  knapsack  problem:  a  dynamic  programming  formulation  and 
an  integer  programming  formulation.  Exact  solution  methods  for  each  formulation 
were  reviewed.  Enumerative  methods  are  typically  applied  to  dynamic  programs, 
while  integer  programs  can  be  relaxed  to  provide  information  about  their  solution. 
Solution  algorithms  that  exploit  the  LP-relaxation  are  the  branch-and-bound  and 
branch-and-cut  algorithms.  Using  the  ideas  associated  with  knapsack  problems,  for¬ 
mal  mathematical  models  were  developed  for  the  payload  selection  problem.  It  was 
assumed  that  a  satellite  constellation  observed  at  fixed  intervals  will  be  maintained 
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or  expanded  by  launching  satellites  individually  at  predetermined  epochs.  At  each 
launch  epoch,  payloads  are  selected  for  the  satellite  bus  to  be  launched  and  assigned 
MMD  specifications.  Each  individual  payload  has  a  utility  associated  with  it  that, 
while  not  satisfying  the  textbook  definition  of  utility,  is  a  function  of  the  payload’s 
importance,  MMD  specification,  and  the  random  number  of  like-type  functional 
payloads.  The  objective  of  the  payload  selection  and  specification  problem  is  to 
maximize  the  total  lifetime  utility  of  the  satellite  constellation.  Four  characteriza¬ 
tions  of  utility  were  discussed:  static  and  deterministic,  dynamic  and  deterministic, 
static  and  stochastic,  and  dynamic  and  stochastic. 

Mathematical  models  for  the  single-satellite  payload  selection  problem  were 
developed  for  each  characterization  of  utility.  It  was  shown  that  each  was  a  relax¬ 
ation  of  a  multi-choice,  multidimensional  knapsack  problem.  Because  dynamic  and 
stochastic  utility  most  closely  describes  the  actual  nature  of  payload  utility,  this 
model  was  extended  to  a  multi-satellite  case  which  was  the  central  focus  of  this 
research.  To  enable  similar  solution  methods  to  those  used  in  knapsack  problems, 
both  dynamic  programming  and  integer  programming  formulations  were  derived  for 
the  multi-satellite  model.  The  dynamic  programming  formulation  was  solved  exactly 
using  a  Matlab  program  to  completely  enumerate  the  state  space  and  prune  infeasi¬ 
ble  states.  The  integer  programming  formulation  was  solved  through  application  of 
Dash  Optimization’s  Xpress  solver  that  employs  the  branch-and-bound  algorithm  in 
concert  with  preprocessing  and  search  heuristics.  In  addition  to  exact  methods,  four 
heuristic  methods  were  introduced.  Each  heuristic  solved  the  payload  selection  and 
specification  problem  as  a  series  of  MCMKPs;  however,  such  an  approach  did  not 
guarantee  optimality.  Both  a  5-norm  and  a  weighted  norm  heuristic  were  developed 
by  extending  the  classical  profit-to-weight  ratio  heuristic  of  the  one-dimensional  KP. 
For  comparison  purposes,  a  greedy  heuristic  that  selects  payloads  based  on  their 
utility  was  introduced.  Finally,  a  simulated  annealing  routine  was  developed  for 
the  payload  selection  problem  using  experimentally  determined  parameters.  Pos- 
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sible  functional  forms  for  payload  utility  functions  and  survival  distributions  were 
developed,  and  the  performance  of  the  exact  and  heuristic  solution  methods  was 
evaluated  on  a  variety  of  both  notional  and  randomly-generated  problem  instances. 
It  was  observed  that  the  Xpress  solver  vastly  outperformed  the  dynamic  program¬ 
ming  enumeration  in  both  solution  time  and  practical  ability  to  attain  an  optimal 
solution.  Of  the  heuristics,  the  solutions  provided  by  simulated  annealing  and  the 
5-norm  heuristic  were  typically  closest  to  optimal.  However,  the  greater  consistency 
of  simulated  annealing  in  providing  near-optimal  solutions  led  to  the  conclusion  that 
it  was  the  best  heuristic  to  use  for  satellite  payload  selection.  Deciding  whether  to 
solve  the  payload  selection  problem  using  an  exact  IP-based  solver  or  a  simulated 
annealing  routine  is  a  multifaceted  decision,  and  tradeoffs  must  be  made  between 
solution  quality  and  complexity. 

Although  a  general  methodology  has  been  developed  for  the  selection  and  spec¬ 
ification  of  satellite  payloads,  there  are  a  number  of  unresolved  issues  and  areas  of 
future  research.  The  power  constraints  of  satellite  buses  were  quantified  in  terms 
of  energy  measurements,  the  lifetime  discharge  of  the  battery.  In  reality,  before  it 
ceases  to  produce  energy,  a  battery’s  true  power ,  the  rate  of  its  energy  delivery,  will 
decline  to  an  insufficient  level  for  the  operation  of  many  payloads.  This  is  a  very 
important  factor  in  the  operability  of  many  satellites,  and  the  model’s  realism  would 
be  greatly  increased  by  incorporating  it.  Additionally,  the  independence  of  payload 
survival  distributions  was  assumed.  This  is  somewhat  restrictive  as  payload  queue¬ 
ing,  in  which  one  payload  prompts  the  operation  of  another,  is  relatively  common. 
Therefore,  the  failure  of  a  payload  may  cause  the  effective  failure  of  another.  Closely 
related  is  extending  the  assumption  of  utility  dependence  to  incorporate  the  utility 
dependence  between  different  types  of  payloads.  A  possible  research  direction  to 
incorporate  dependence  in  survival  distributions  and  utility  functions  is  attempting 
to  aggregate  the  dependent  payloads  into  a  single  payload  whose  survival  distribu¬ 
tion  and  utility  function  characterize  the  superposition  of  the  aggregated  payloads. 
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Finally,  neither  thermal  or  geometric  constraints  were  considered.  In  practice  it  is 
critical  to  dissipate  the  heat  generated  by  payloads.  Also,  it  must  be  ensured  that  the 
payload’s  geometries  allow  them  to  fit  properly  on  the  bus.  The  ability  to  extend  the 
payload  selection  and  specification  model  to  incorporate  these  items  would  greatly 
enhance  both  its  usefulness  and  effectiveness  in  the  process  of  acquiring,  launching 
and  operating  satellite  systems. 
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Appendix  A.  DP  Enumeration  Code 

%7.%7o%y.%7.%yo%%7//o7oyo%%7.%7o%y.%7.%yo%%%7.%yo%%%7.7o%0/.%7.7.yo7.0/.%7.7o0/o0/.%%7o7o7.0/.%7.%yo%0/.%7.7o7o%%7.7o7oyo°/.% 

#/o  AUTHOR:  Capt  John  Flory 
l  AFIT/ENS/ G0R-06M 

°/0  March  2006 

°/0  This  program  performs  a  complete  enumeration  of  the  multi-satellite 
°/0  payload  selection  and  specification  problem.  Pruning  is  used  to 
l  eliminate  infeasible  states. 

y.mmmmmrammmraramrara%mmmmmmmramnn 


clear  all; 

°/0  Begin  clock 
tic ; 

°/0  Set  alpha  and  beta,  the  tuning  parameters  for  the  utility  decay  function 
°/0  and  survival  distribution,  respectively 
alpha  =  log( . 5) ; 
beta  =  log( . 9) ; 

gamma  =  . 5 ; 

°/0  Number  of  satellites,  payload  types  and  MMD  specifications  per  type 
N_Sat  =  2; 

N_Type  =  5; 

N_Spec  =  3; 

°/0  Time  horizion  of  problem 
Epochs  =  30; 

l  Set  launch  time  periods 
N1  =  [2  5  7  9  10  11  13  15] ; 


%  Input  Payload  Data  —  [Importance,MMD,Power,Cost,Weight, Volume] 


PD(1 , : 

, 1)= [10 

3 

500 

425 

450 

15] 

PD (2,  : 

, 1)= [10 

6 

500 

460 

475 

17] 

PD (3,  : 

, 1)= [10 

10 

500 

500 

500 

20] 

PD(1 , : 

,2)  =  [8.5 

3 

475 

375 

400 

16]  ; 

PD (2,  : 

,2)  =  [8.5 

6 

475 

405 

415 

18]  ; 

PD (3,  : 

,  2)  =  [8 . 5 

10 

475 

430 

420 

19. 

PD(1 , : 

,3)  =  [7.5 

3 

425 

410 

430 

10]  ; 

42 

PD(2, 

,  3)  =  [7 .5 

6 

425 

460 

480 

13]  ; 

43 

PD  (3, 

,3)  =  [7 .5 

10 

425 

480 

495 

14]; 

44 

PD(1 , 

,4)  =  [7 

3 

260 

300 

230 

10]; 

45 

PD  (2, 

,4)  =  [7 

6 

260 

350 

280 

13]  ; 

46 

PD  (3, 

,4)  =  [7 

10 

260 

370 

300 

14]; 

47 

PD(1 , 

,5)  =  [6 

3 

225 

370 

380 

13]  ; 

48 

PD  (2, 

,5)  =  [6 

6 

225 

400 

390 

15.5]  ; 

49 

PD  (3, 

,5)  =  [6 

10 

225 

410 

395 

17.5]  ; 

50 

PD(1 , 

,6)  =  [5.5 

3 

300 

280 

240 

8]; 

51 

PD  (2, 

,6)  =  [5.5 

6 

300 

320 

290 

9]; 

52 

PD  (3, 

,  6)  =  [5 . 5 

10 

300 

380 

310 

12]  ; 

53 

PD(1 , 

,7)  =  [5 

3 

275 

150 

280 

7]  ; 

54 

PD  (2, 

,7)  =  [5 

6 

275 

190 

350 

9.5]  ; 

55 

PD  (3, 

,7)  =  [5 

10 

275 

240 

410 

14]; 

56 

PD(1 , 

,  8)  =  [3 

3 

175 

270 

225 

4]; 

57 

PD  (2, 

,8)  =  [3 

6 

175 

310 

260 

5.5]  ; 

58 

PD  (3, 

,8)  =  [3 

10 

175 

335 

300 

8]; 

59 

60  7,  Create  array  of  the  initial  number  of  each  payload  type  in  constellation 

61  7,  at  time  0 

62  Init_Cons  =  zeros (N_Type ,N_Spec) ; 

63 

64  7«  Input  the  number  of  each  specific  payload  type/nonzero  MMD  combination  in 

65  7.  constellation  at  time  0 

66  Init_Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;0  0  0] ; 

67 

68  7*  Create  array  for  expected  numbers  of  payload  types  at  each  time  period 

69  ENum_0rig  =  zeros (N_Type , Epochs) ; 

70 

71  7o  Given  the  initial  numbers  of  each  payload  type/nonzero  MMD  combination 

72  7o  initially  in  constellation,  calculate  their  expected  numbers  at  future 

73  7*  time  periods. 

74  for  i  =  1: Epochs; 

75  for  j  =  l:N_Type; 

76  Num  =  0; 

77  for  k  =  l:N_Spec 

78  if  i  <=  PD(k,2,j)+l; 

79  Num  =  Num  +  Init_Cons ( j ,k)  *  Death_Exp(beta,PD(k,2, j) ,i-l) ; 

80  end; 

81  end; 

82  ENum_0rig( j , i)  =  Num; 

83  end ; 

84  end ; 
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85 


86 

87 
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7,  Set  number  of  states  for  stage  1 
N_States  =  N_Spec+l; 

N_States 


89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 
127 


7*  Initialize  utility  vector  for  all  states  in  stage  1 
Utility  =  zeros (N_States , 1) ; 

7*  Initialize  utility  vector  for  all  states  in  previous  stage 
Utility_Last  =  zeros (N_States , 1) ; 

7,  Initialize  vector  that  stores  payload  contents  of  stage  1 
Contents  =  zeros (N_States , 1) ; 

7,  Initialize  vector  of  expected  number  of  payload  type  being  specified 
7,  in  stage  1 

ENum  =  zeros (N_States , Epochs) ; 

7,  Initialize  vector  of  state  designators  as  feasible  or  infeasible 
Feasible  =  zeros (N_States , 1) ; 

7,  Initialize  vectors  of  satellites ’  remaining  power,  costs,  weight,  and 
7,  volume  resources  for  all  states  in  stage  1 
Power  =  zeros (N_States ,N_Sat) ; 

Cost  =  zeros (N_States ,N_Sat) ; 

Weight  =  zeros (N_States ,N_Sat) ; 

Volume  =  zeros (N_States ,N_Sat) ; 

7,  Set  power,  cost,  weight,  and  volume  capacities  of  all  satellite  buses 
for  i  =  l:N_States; 
for  j  =  l:N_Sat; 

Power(i,j)  =  10500; 

Cost(i,j)  =  2500; 

Weight (i,j)  =  2500; 

Volume (i,j)  =  100; 

end; 

end; 

7,  Create  identity  matrix  to  construct  elementary  vectors 
Ident  =  eye(N_Sat ,N_Sat) ; 

7*  Begin  stage  1 
for  i  =  l:N_Type; 
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for  j  =  l:N_Sat; 

e_i  =  Ident ( j , : ) ; 

/  Only  do  for  stage  1 
if  i==l  &  j==l ; 

7,  Iterate  through  all  stage  1  states 
for  k  =  l:N_States; 

l  Load  expected  number  of  payload  type  being  added 
7,  stage  1 

ENum(k,:)  =  ENum_0rig(i , : ) ; 

7#  Update  contents  vector  to  reflect  the  addition  of  this 
7o  payload  type  and  its  given  MMD  specification 
Contents  (k,  1)  =  mod(k,N_Spec+l)  *10~  (N_Sat* (i-1) +  (  j-1) )  ; 

7#  If  specifying  a  nonzero  MMD 
if  mod(k,N_Spec+l)  ~=  0; 

MMD  =  PD(k,2,i); 

psi  =  PD(mod(k,N_Spec+l) , 1 , i) ; 

7.  Iterate  over  all  epochs  payload  is  in  service 
for  1  =  Nl(j)+1  :  (N1 (j )+l)+MMD 
n  =  1-N1 ( j ) -1 ; 

7.  Update  number  of  expected  payloads  at  each 
7.  time 

ENum(k,l)=  . . . 

ENum_0rig(i , l)+l*Death_Exp(beta,MMD,n) ; 

7.  Get  total  number  of  like-type  payloads  in 
7#  the  constellation 
N  =  max(ENum(k,l) , 1) ; 

l  Compute  the  power  scaling  factor 
a  =  PD(mod(k,N_Spec+l) ,3, i)/ . . . 

Util_Exp(psi , 1 , gamma, alpha, MMD, 0) ; 

7.  Compute  the  added  utility  of  a  time 
i  period 

z  =  Ut il_Exp (psi ,N, gamma, alpha, MMD, n) .. . 
*Death_Exp(beta,MMD,n) ; 
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l  Update  total  payload  utility 
/  and  bus  power  consumption 
Utility(k,l)  =  Utility(k,l)  +  z; 

Power (k,:)  =  Power (k,:)  -  z  *  a  *  e_i; 

end; 

end; 

i  If  payload  has  nonzero  MMD  specification  subtract 
°/0  power,  cost,  weight,  and  volume  resources  from  bus 
°/0  capacity 

if  mod(k,N_Spec+l)  ~=  0; 

Cost(k,:)  =  Cost(k, : )-PD(mod(k,N_Spec+l) ,4,i)*e_i; 
Weight(k,:)  =  Weight(k, : )-PD(mod(k,N_Spec+l) ,5,i)*e_i; 
Volume(k,:)  =  Volume(k, : )-PD(mod(k,N_Spec+l) ,6, i)*e_i ; 

end; 

end; 

7,  Now  check  feasibility  of  states  in  stage  and  count  number  of 
%  feasible  states 
N_Feasible  =  N_States; 

for  k  =  l:N_States 

if  Power (k,j)<0  I  Cost(k,j)<0  I  Weight (k,j)<0  I  ... 

Volume (k,j)<0 
Feasible (k, 1)  =  0; 

N_Feasible  =  N_Feasible-l ; 

else 

Feasible (k, 1)  =  1; 

end; 

end; 

7«  Now  eliminate  infeasible  states 
t  =  1; 

for  k  =  l:N_States 

if  Feasible(k, 1)  ==  1; 

Utility(t,l)  =  Utility (k, 1) ; 

Contents (t , 1)  =  Contents (k, 1) ; 

Power (t,:)  =  Power(k,:); 

Cost(t,:)  =  Cost(k,:); 

Weight (t,:)  =  Weight(k,:); 

Volume (t,:)  =  Volume(k,:); 

ENum(t,:)  =  ENum(k,:); 
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214 


t  =  t+1; 


215 

216 
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235 
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end; 

end; 

7,  Transfer  information  about  current  states  into  storage  to 

%  transfer  to  stage  2 

N_States  =  N_Feasible  *  (N_Spec+l) ; 

Utility_Last  =  Utility; 

Contents_Last  =  Contents; 

ENum_Last  =  ENum; 

Power_Last  =  Power; 

Cost_Last  =  Cost; 

Weight_Last  =  Weight; 

Volume_Last  =  Volume; 

N_States 


else 

°i  Iterate  through  all  states  in  stage  and  load  state  values 
7«  from  the  previous  stage 
for  k  =  l:N_States; 

Utility (k,l)  =  Utility_Last (ceil (k/ (N_Spec+l) ) , 1) ; 
Contents (k, 1)  =  Contents_Last (ceil(k/(N_Spec+l) ) , 1) . . . 

+  mod(k,N_Spec+l)  *  10" (N_Sat* (i-l)+( j-1) ) ; 

Power (k,:)  =  Power_Last (ceil (k/ (N_Spec+l) ) , : ) ; 

Cost(k,:)  =  Cost_Last(ceil(k/(N_Spec+l)) , :) ; 

Weight (k,:)  =  Weight_Last (ceil (k/ (N_Spec+l) ) , : ) ; 

Volume (k,:)  =  Volume_Last (ceil (k/ (N_Spec+l) ) , : ) ; 

7»  If  you  start  specifying  next  payload  type,  load  new 
l  expected  numbers;  otherwise,  continue  with  the  old. 

if  j  '=  l; 

ENum(k,:)  =  ENum_Last (ceil (k/ (N_Spec+l) ) , : ) ; 
else ; 

ENum(k,:)  =  ENum_0rig(i , : ) ; 

end; 


end; 

#/«  Iterate  through  all  states  in  stage 
for  k  =  l:N_States; 

l  If  payload  has  nonzero  MMD 
if  mod(k,N_Spec+l)  ~=  0; 
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7,  Load  payload  data  values 
MMD  =  PD(mod(k,N_Spec+l) ,2,i) ; 
psi  =  PD(mod(k,N_Spec+l) , 1 , i) ; 

7*  Compute  no.  time  periods  payload  is  in  service 
for  1  =  Nl(j)  +  1  :  (Nl(j)+1)  +  MMD 
n  =  1  -  Nl(j)  -  1; 

7o  Update  expected  number  of  functional  payloads 
ENum(k,l)  =  ENum(k,l)  +  Death_Exp (beta, MMD, n) ; 

N  =  max(ENum(k,l) , 1) ; 

7.  Compute  power  scaling  factor 
a  =  PD(mod(k,N_Spec+l) ,3, i)/ .  .  . 

Util_Exp(psi , 1 , gamma, alpha, MMD, 0) ; 

7o  Compute  stage  incremental  utility 
z  =  Ut il_Exp (psi ,N, gamma, alpha, MMD, n) .. . 

*  Death_Exp(beta,MMD,n) ; 

7.  Update  total  payload  utility  and 
7o  bus  power  consumption 
Utility(k,l)  =  Utility(k,l)  +  z; 

Power (k,:)  =  Power (k,:)  -  a  *  z  *  e_i; 

end; 

end; 

7.  Update  total  bus  power,  weight,  and  volume  consumption 
if  mod(k,N_Spec+l)  ~=  0; 

Cost(k,:)  =  Cost(k,:)  -PD(mod(k,N_Spec+l) ,4,i)*e_i; 
Weight(k,:)  =  Weight(k,:)  -PD(mod(k,N_Spec+l) ,5,i)*e_i 
Volume(k,:)  =  Volume(k,:)  -PD(mod(k,N_Spec+l) ,6,i)*e_i 

end; 

end; 

7,  Now  check  feasibility  of  states  in  stage  and  count  number  of 
7.  feasible  states 
N_Feasible  =  N_States; 

for  k  =  l:N_States 

if  Power (k,j)<0  I  Cost(k,j)<0  I  Weight (k,j)<0. . . 

I  Volume (k,j)<0 
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Feasible (k, 1)=0 ; 

N_Feasible  =  N_Feasible-l ; 

else 

Feasible (k, 1)  =  1; 

end; 

end; 

7,  Now  eliminate  infeasible  states 
t  =  1; 

for  k  =  l:N_States 

if  Feasible(k, 1)  ==  1; 

Utility(t,l)  =  Utility (k, 1) ; 

Contents (t , 1)  =  Contents (k, 1) ; 

Power (t,:)  =  Power(k,:); 

Cost(t,:)  =  Cost(k,:); 

Weight (t,:)  =  Weight(k,:); 

Volume (t,:)  =  Volume(k,:); 

ENum(t,:)  =  ENum(k,:); 
t  =  t+1; 

end; 

end; 

if(i~=N_Type  I  j ~=N_Sat) ; 

l  Set  number  of  states  for  next  stage 
N_States  =  N_Feasible  *  (N_Spec+l) ; 

else 

N_States  =  N_Feasible; 

end; 

7«  Store  all  stage  values  for  use  in  next  stage 
Utility_Last  =  Utility; 

Contents_Last  =  Contents; 

ENum_Last  =  ENum; 

Power_Last  =  Power; 

Cost_Last  =  Cost; 

Weight_Last  =  Weight; 

Volume_Last  =  Volume; 

N_States 

end; 

end; 

end; 

disp(’ENDO  ; 
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7o  Find  end  state  with  greatest  total  utility 
Best  =  0; 

Best_Index  =  1; 
for  k  =  l:N_States; 

if  Utility (k,l)  >  Best 
Best  =  Utility (k, 1) ; 

Best_Index  =  k; 

end; 

end; 

Sat_Cont  =  zeros (N_Sat ,N_Type) ; 

Best_Content  =  Contents (Best_Index , 1) ; 

7,  Translate  the  contents  of  the  best  state  into  actual  payload 
7,  specifications 
for  i  =  l:N_Type; 

for  j  =  l:N_Sat; 

z  =  floor (Best_Content/10'‘ (  (N_Sat* (N_Type-l)  +.  .  . 

(N_Sat-l) )  -  N_Sat*(i-l)-(j-l))) ; 

Sat_Cont (N_Sat+l-j ,N_Type+l-i)  =  z; 

Best_Content  =  Best_Content 

z*10~((N_Sat*(N_Type-l)+(N_Sat-l))-N_Sat*(i-l)-(j-l)) ; 

end; 

end; 

7,  Stop  clock 
toe ; 

time  =  toe; 
time 

7.  Display  solution,  optimal  utility  value  and  remaining  power,  cost, 

7,  weight  and  volume  on  all  buses 

Sat_Cont 

Utility (Best _Index) 

Power (Best_Index, :) 

Cost (Best_Index , : ) 

Weight (Best_Index , : ) 

Volume (Best_Index , : ) 
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Appendix  B.  IP  Generation  Code 


y.7.y;/.7.y:/:/.y:/:/.y.y:/;/.y;/;/.7.y//.7.y:/:/.y:/:/.y.y:/.y.y:/;/.y:/;/.7.y:/.7.y//:/.y//:/.y.y:/;/.y//;/.%y;/.7.y:/:/.y.,/.7.y.y:/.y.y:/.y.y:/:/. 

#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 

7,  March  2006 

7,  This  program  generates  the  objective  functions  and  constraint  matrix  for 
7*  the  IP  formulation  of  the  payload  selection  and  specification  problem. 
7o7o707o7.7»7o7«7«7o7o7o7.7«7o7.7o7o7o7o7o7o7»7.7«7o7o7o7o7.7o7o7o7o7o7.707o7.7#7o7«7o7o7o7o7.7«7o7.7o7o7o7o7o7.7#7.7«7o7.7«7o7.7o7o7.7o7o7o7o7o7o7o7o 

7.  Set  survival  and  utility  decay  parameters 
Death_Coef  =  log(.9); 

Util_Coef  =  log(.5); 

7,  Set  utility  dependence  parameter 
Type_Dep  =  .5; 

7,  Set  number  of  satellites,  payload  types,  and  nonzero  MMD  specifications 
N_Sat  =  2; 

N_Type  =  5; 

N_Spec  =  3; 

7,  Time  horizon  of  problem 
Epochs  =  30; 

7,  Set  launch  time  periods 
N1  =  [2  5  7  9  10  11  13  15] ; 


27 

7,  Input  Payload 

Data  — 

[Importance , MMD , Power , Cost 

Weight , 

28 

PD(1 , : 

II 

o 

3 

500 

425 

450 

15]  ; 

29 

PD (2,  : 

II 

o 

6 

500 

460 

475 

17]; 

30 

PD (3, : 

II 

o 

10 

500 

500 

500 

20]  ; 

31 

PD(1, : 

to 

II 

00 

01 

3 

475 

375 

400 

16]  ; 

32 

PD (2,  : 

to 

II 

00 

cn 

6 

475 

405 

415 

18]  ; 

33 

PD (3, : 

to 

II 

00 

CJ1 

10 

475 

430 

420 

19.5 

34 

PD(1, : 

CO 

II 

CJ1 

3 

425 

410 

430 

10]  ; 

35 

PD (2,  : 

CO 

II 

CJ1 

6 

425 

460 

480 

13]  ; 

36 

PD (3, : 

CO 

II 

01 

10 

425 

480 

495 

14]; 

37 

PD(1 , : 

n 

3 

260 

300 

230 

10]  ; 

38 

PD (2,  : 

n 

6 

260 

350 

280 

13]  ; 

39 

PD (3, : 

n 

10 

260 

370 

300 

14]; 

40 

PD(1, : 

,5)  =  [6 

3 

225 

370 

380 

13]  ; 

PD  (2, 


, 5)  =  [6 


15.5] ; 


41 


6 


225 


400 


390 


42 

PD(3, 

:  ,5)  =  [6 

10 

225 

410 

395 

17.5]  ; 

43 

PD(1 , 

:  ,6)  =  [5.5 

3 

300 

280 

240 

8]; 

44 

PD  (2, 

:  ,6)  =  [5.5 

6 

300 

320 

290 

9]; 

45 

PD  (3, 

:  ,6)  =  [5.5 

10 

300 

380 

310 

12]  ; 

46 

PD(1 , 

:  ,7)  =  [5 

3 

275 

150 

280 

7]; 

47 

PD  (2, 

:  ,7)  =  [5 

6 

275 

190 

350 

9.5]  ; 

48 

PD  (3, 

:  ,7)  =  [5 

10 

275 

240 

410 

14]; 

49 

PD(1, 

:,8)  =  [3 

3 

175 

270 

225 

4]; 

50 

PD  (2, 

:,8)  =  [3 

6 

175 

310 

260 

5.5]  ; 

51 

PD  (3, 

:,8)  =  [3 

10 

175 

335 

300 

8]; 

52 

53 

"/,  Set 

satellite  bus 

capacities 

54 

Power 

_Limit  =  10500 

55  Cost_Limit  =  2500; 

56  Volume_Limit  =  2500; 

57  Weight_Limit  =  100; 

58  ENum_0rig  =  zeros (N_Type , Epochs) ; 

59 

60  7,  Set  numbers  of  payload  types/MMD  specifications  initially  in  the 

61  /  constellation 

62  Init_Cons  =  zeros (N_Type ,N_Spec) ; 

63  Init_Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0] ; 

64 

65  7,  Calculate  expected  number  of  remaining  payloads  at  each  time  period 

66  for  i  =  1: Epochs; 

67  for  j  =  l:N_Type; 

68  Num=0 ; 

69  for  k  =  l:N_Spec 

70  if  i  <=  PD(k,2,j)  +  1; 

71  Num  =  Num  +  Init_Cons ( j ,k)  *  ... 

72  Death_Exp(Death_Coef ,PD(k,2, j) ,i-l) ; 

73  end ; 

74  end ; 

75  ENum_0rig(j , i)  =  Num; 

76  end ; 

77  end ; 

78 

79  7«  Compute  total  number  of  variables  and  number  of  variables  associated  with 

80  7o  each  payload  type 

81  tot_var  =  N_Type  *  (N_Spec  +  l)~N_Sat; 

82  type_var  =  (N_Spec  +  l)~N_Sat; 

83 

84  7«  Create  vector  of  utility  objective  coefficients,  and  power,  cost,  weight, 
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86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 
127 


7,  and  volume  coefficients  for  the  constraints 
Utility  =  zeros (1 ,tot_var) ; 

Power  =  zeros (tot_var,  N_Sat) ; 

Cost  =  zeros (tot_var,  N_Sat) ; 

Volume  =  zeros (tot_var,  N_Sat) ; 

Weight  =  zeros (tot_var,  N_Sat) ; 

7,  Create  array  of  all  possible  MMD  specifications  over  the  satellites  for 

7,  each  payload  type 

Spec  =  zeros (tot_var,N_Sat) ; 

for  i  =  1  :  N_Type; 

for  j  =  1  :  type_var; 

7,  Translate  variable  index  into  a  set  of  MMD  specifications 
comb  =  dec2base( j-1 ,4) ; 
dec_comb  =  str2num(comb) ; 
spec  =  zeros(l ,N_Sat) ; 

for  k  =  1  :  N_Sat; 

ENum(k,:)  =  ENum_Orig(i , : ) ; 

end; 

7,  Generate  MMD  specifications 
for  k  =  1  :  N_Sat; 

spec(l, k)  =  floor (dec_comb/10~ (N_Sat-k) ) ; 
dec_comb  =  dec_comb  -  spec(l, k)  *  10" (N_Sat-k) ; 

end; 

Spec(type_var  *  ( i— 1 )  +  j,:)  =  spec; 
u  =  0;  ul  =  0; 


7,  Calculate  expected  number  of  payloads 
for  k  =  1  :  N_Sat; 

posn  =  N_Sat  -  (k— 1 ) ; 
if  speed, posn)  ~=  0; 

MMD  =  PD(spec(l ,posn) ,2,i) ; 
for  1  =  1  :  MMD  +  1; 

ENum(l ,Nl(posn)+l)  =  ENum(l ,Nl(posn)+l)  +  ... 
Death_Exp(Death_Coef , MMD, 1-1) ; 

end; 

end; 

end; 
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°/0  Calculate  utility,  power,  cost,  weight,  volume  requirements  for 

129 

7,  each  set  of  MMD  specifications 

130 

for  k  =  1  :  N_Sat; 

131 

p  =  0 ;  c  =  0 ;  w  =  0 ;  v  =  0 ;  ut  =  0; 

132 

posn  =  N_Sat  -  (k-1) ; 

133 

if  spec(l, posn)  ~=  0; 

134 

psi  =  PD(spec(l ,posn) , 1 ,i) ; 

135 

MMD  =  PD(spec(l ,posn) ,2,i) ; 

136 

a_l  =  PD(spec(l ,posn) ,3,i)  /  ... 

137 

Util_Exp(psi , 1 ,Type_Dep,Util_Coef ,MMD,0) ; 

138 

c  =  c  +  PD(spec (1 ,posn) ,4, i) ; 

139 

w  =  w  +  PD(spec (1 ,posn) ,5,i) ; 

140 

v  =  v  +  PD (spec d , posn) ,6,i) ; 

141 

142 

7,  Compute  total  utility  and  power  requirements 

143 

for  1  =  1  :  MMD  +  1; 

144 

N  =  max(  ENum(l ,N1 (posn)+l) , 1) ; 

145 

ul  =  Util_Exp(  psi ,N,Type_Dep,Util_Coef , MMD, 1-1  )  *  ... 

146 

Death_Exp(  Death_Coef ,MMD , 1-1  ); 

147 

u  =  u  +  ul; 

148 

p  =  p  +  ul  *  a_l; 

149 

ENum(l ,Nl(posn)+l)  =  ENum(l ,Nl(posn)+l)  -  ... 

150 

Death_Exp(Death_Coef , MMD, 1-1) ; 

151 

ut  =  ut+ul ; 

152 

end; 

153 

154 

l  Set  all  coefficients  with  computed  values 

155 

end; 

156 

UtilityK  type_var*(i-l)+j  ,posn  )=ut; 

157 

Power(  type_var*(i-l)+j ,  posn  )  =  p; 

158 

Cost(  type_var*(i-l)+j ,  posn  )  =  c; 

159 

Volume(  type_var* (i-l)+j ,  posn  )  =  v; 

160 

Weight (  type_var* (i-l)+j ,  posn  )  =  w; 

161 

end; 

162 

Utility(l,  type_var*(i-l)+j  )  =  u; 

163 

end; 

164 

end; 

165 

166 

/  Create  vector  of  objective  coefficients 

167 

0bj_Coefs  =  Utility; 

168 

169 

7,  Compute  number  of  constraints 

170 

num_constraint  =  N_Type  +  4  *  N_Sat; 
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l  Create  constraint  matrix  and  RHS  constraints 
Cons_Mat  =  zeros (num_constraint ,  tot_var) ; 

RHS  =  zeros (num_constraint , 1) ; 

°/0  Create  special  ordered  set  constraints 
for  i  =  1  :  N_Type; 

for  j  =  1  :  type_var; 

Cons_Mat(i,  type_var* ( i— 1 )  +  j)  =1; 

RHS(i)  =  1; 

end; 

end; 

#/o  Create  power,  cost,  weight,  and  volume  constraints 
for  i  =  1  :  N_Type; 

for  j  =  1  :  type_var; 
for  k  =  1  :  N_Sat; 

cons  =  (N_Type  +  1)  +  4  *  (k-1) ; 
indx  =  type_var  *  ( i— 1 )  +  j; 

Cons_Mat (cons , indx)  =  Power (indx, N_Sat- (k-1) ) ; 
Cons_Mat (cons+1 , indx)  =  Cost (indx, N_Sat- (k-1) ) ; 
Cons_Mat(cons+2,indx)  =  Weight (indx, N_Sat- (k-1) ) ; 
Cons_Mat(cons+3,indx)  =  Volume (indx, N_Sat- (k-1) ) ; 
RHS (cons)  =  Power_Limit; 

RHS (cons+1)  =  Cost_Limit; 

RHS(cons+2)  =  Volume_Limit ; 

RHS(cons+3)  =  Weight_Limit ; 

end; 

end; 

end; 

Cons_Matl  =  zeros (num_ constraint ,  tot_var+l) ; 

for  i  =  1  :  num_constraint ; 

Cons_Matl (i , : )  =  [Cons_Mat (i , : ) , -1] ; 

end; 

7,  Create  array  with  just  objective  and  constraint  coefficients 
Systeml  =  [Utility,  -1 ; Cons_Matl] ; 

Systeml  =  Systeml’; 

7.  Output  array  to  .csv  file 
csvwrite(’ Systeml.txt’ , Systeml) ; 
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Appendix  C.  Simulated  Annealing  Code 


y.7.y:/.7.y:/:/.y:/:/.y.y:/;/.y:/;/.7.y:/.7.y:/:/.y:/:/.y.y:/.y.y:/:/.y:/;/.7.y:/.7.y:/:/.y//:/.y.y:/.y.y:/.y.7.,/:/.7.,/:/.7.,/:/.7.y.y.7.y.y:/.y.y:/:/. 


#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 


7,  March  2006 

7,  This  program  applies  a  simulated  annealing  routine  to  the  payload 
7,  specification  problem.  Simulated  annealing  is  used  to  maximize  the 
7,  utility  of  each  satellite  in  succession  by  solving  the  multi-choice, 

7,  multidimensional  knapsack  problem  associated  with  each  satellite. 

7,  Multiple  replications  of  this  proceedure  are  applied  to  the  constellation 
7,  and  the  replication  with  the  greatest  overall  utility  is  chosen.  A 
7,  geometric-cooling  schedule  is  used. 


7o7.7o7o7.7o7o7.7o7o7.7o7«7o7o7o7o7o7.7o7o7.7o7o7.7o7o7o7o7«7.7o7o7o7o7o7o7o7.7o7o7.7o7o7.7o7o7.7o7o7o7o7o707o7.7o7o7.707o7.7o7o7.7o7o7.7o7o7o7o7o7o7o 


clear  all; 
tic ; 

7,  Set  number  of  simulated  annealing  replications 
N_replicate  =  30; 

7,  Set  survival  and  utility  function  tuning  parameters 
Death_Coef  =  log(.9); 

Util_Coef  =  log(.5); 

7,  Set  utility  dependence  parameter 
Type_Dep  =  .5; 

7,  Set  number  of  satellites,  payload  types,  and  nonzero  MMD  specifications 
7,  available  to  each  payload  type 
N_Sat  =  7; 

N_Type  =  8; 

N_Spec  =  3; 

7,  Set  time  horizon  of  problem 
Epochs  =  30; 

7,  Set  satellite  launch  time  periods 
N1  =  [2  5  7  9  10  11  13  15] ; 

7,  Input  Payload  Data  —  [Importance,MMD,Power,Cost,Weight, Volume] 

PD(1 , : , 1)  =  [10  3  500  425  450  15]; 

PD(2 , : , 1)  =  [10  6  500  460  475  17]; 
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42 

PD(3, : 

II 

o 

10 

500 

500 

500 

20]  ; 

43 

PD(1 , : 

to 

II 

00 

cn 

3 

475 

375 

400 

16]  ; 

44 

PD (2,  : 

to 

II 

00 

CJ1 

6 

475 

405 

415 

18]  ; 

45 

PD (3, : 

to 

II 

00 

CJ1 

10 

475 

430 

420 

19.5]  ; 

46 

PD(1 , : 

CO 

II 

01 

3 

425 

410 

430 

10]; 

47 

PD (2,  : 

,3)  =  [7.5 

6 

425 

460 

480 

13]  ; 

48 

PD (3, : 

,3)  =  [7.5 

10 

425 

480 

495 

14]; 

49 

PD(1 , : 

,4)  =  [7 

3 

260 

300 

230 

10]; 

50 

PD (2,  : 

,4)  =  [7 

6 

260 

350 

280 

13]  ; 

51 

PD (3, : 

,4)  =  [7 

10 

260 

370 

300 

14]; 

52 

PD(1 , : 

,5)  =  [6 

3 

225 

370 

380 

13]  ; 

53 

PD (2,  : 

,5)  =  [6 

6 

225 

400 

390 

15.5]  ; 

54 

PD (3, : 

,  5)  =  [6 

10 

225 

410 

395 

17.5]  ; 

55 

PD(1, : 

,6)  =  [5.5 

3 

300 

280 

240 

8]; 

56 

PD (2,  : 

,6)  =  [5.5 

6 

300 

320 

290 

9]; 

57 

PD (3, : 

,6)  =  [5.5 

10 

300 

380 

310 

12]  ; 

58 

PD(1 , : 

,7)  =  [5 

3 

275 

150 

280 

7]  ; 

59 

PD (2,  : 

,7)  =  [5 

6 

275 

190 

350 

9.5]  ; 

60 

PD (3, : 

,7)  =  [5 

10 

275 

240 

410 

14]; 

61 

PD(1 , : 

,8)  =  [3 

3 

175 

270 

225 

4]; 

62 

PD (2,  : 

,8)  =  [3 

6 

175 

310 

260 

5.5]  ; 

63 

PD (3, : 

,8)  =  [3 

10 

175 

335 

300 

8]; 

64 

65  7,  Input  satellite  bus  capacities 

66  for  i  =  l:N_Sat; 

67  Power(i)  =  10500; 

68  Cost(i)  =  2500; 

69  Weight (i)  =  2500; 

70  Volume (i)  =  100; 

71  end ; 

72 

73  7*  Create  empty  arrays  to  store  the  solution  and  total  utility  of  each 

74  7«  replication 

75  rep_soln  =  zeros (N_Sat ,N_Type ,N_replicate) ; 

76  rep_utility  =  zeros(l,N_replicate) ; 

77 

78  7o  For  all  replications 

79  for  replicate  =  1  :  N_replicate; 

80  replicate 

81  ENum_0rig  =  zeros (N_Type , Epochs) ; 

82 

83  7«  Input  numbers  of  each  payload  type/MMD  specification  initially  in 

84  7o  the  constellation 


C-2 


85 


86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 
127 


Init_Cons  =  zeros (N_Type ,N_Spec) ; 

Init.Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0] ; 

°/«  Calculate  expected  number  of  remaining  payloads  from  those 
°/0  originally  in  the  constellation 
for  i  =  1: Epochs; 

for  j  =  l:N_Type; 

Num=0 ; 

for  k  =  l:N_Spec 

if  i  <=  PD(k,2,j)  +  1; 

Num  =  Num  +  Init_Cons ( j ,k)  *  Death_Exp(Death_Coef ,PD(k, 2 , j ) , i-1) 

end; 

end; 

ENum(j,i)  =  Num; 

end; 

end; 

7«  Create  arrays  to  store  total  utility  of  each  satellite  and  the 
°/0  specifications  of  its  payloads 
Utility  =  zeros(l ,N_Sat) ; 

Soln  =  zeros (N_Sat ,N_Type) ; 

for  i  =  l:N_Sat; 

7,  Create  vectors  to  store  the  utility,  power,  cost,  weight,  and 

/  volumes  of  all  payload/specif ication  combinations 

utility  =  zeros(l,  N_Type  *  N_Spec) ; 

power  =  zeros (1,  N_Type  *  N_Spec) ; 

cost  =  zeros(l,  N_Type  *  N_Spec) ; 

weight  =  zeros (1,  N_Type  *  N_Spec) ; 

volume  =  zeros (1,  N_Type  *  N_Spec) ; 

for  j  =  1  :  N_Type; 

for  k  =  1  :  N_Spec; 

7#  Compute  cost,  weight,  and  volume  of  payload/specif  ication 
7.  combination 

cost(l,  N_Spec  *  ( j  —  1 )  +  k  )  =  PD(k,4,j); 

weight (1,  N_Spec  *  (j-1)  +  k  )  =  PD(k,5,j); 

volume (1,  N_Spec  *  (j-1)  +  k  )  =  PD(k,6,j); 

7#  Load  importance  and  MMD  of  combination 
psi  =  PD(k, 1 , j) ; 

MMD  =  PD(k,2, j) ; 
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•/.  Initialize  utility  and  power  summation  values  to  zero 
u  =  0;  ul  =  0; 

P  =  0; 

/  Compute  power  scaling  factor 

a_l  =  PD(k,3, j)  /  Util_Exp(psi , 1 ,Type_Dep,Util_Coef ,MMD,0) ; 

/  Compute  total  utility  and  power  consumption  of 
/  payload/specif ication  combination 
for  1  =  1  :  MMD  +  1; 


7#  Compute  expected  number  of  payloads  if  payload  is 
7,  included 

N  =  max(  ENum( j ,N1 (i)+l)  +  ... 

Death_Exp(Death_Coef ,  MMD,  1-1),  1); 
ul  =  Util_Exp(  psi ,N,Type_Dep,Util_Coef , MMD, 1-1  )  ... 

*  Death_Exp(  Death_Coef , MMD, 1-1  ); 
u  =  u  +  ul; 
p  =  p  +  ul  *  a_l; 

end; 

7#  Store  total  utility  and  power  consumption  of  combination 
utility (1,  N_Spec  *  ( j  —  1 )  +  k  )  =  u; 
power (1,  N_Spec  *  ( j  —  1 )  +  k  )  =  p; 

end; 

7,  Create  utility/power/cost/weight/volume  master  array 
UPCWV  =  [utility’ , power’ , cost ’ , weight’ , volume ’] ; 


end; 

7,  Randomly  generate  a  starting  slut  ion 
Vc  =  zeros (1 ,N_Type*N_Spec) ; 

for  j  =  1  :  N_Type; 

7,  Decide  whether  to  include  (1)  or  exclude  (0)  payload  type 
inc  =  round (rand(l , 1) ) ; 

7,  If  including,  determine  random  specification  for  payload 
if  inc  ==  1; 

inc_spec  =  round (N_Spec  *  rand (1,1)  -  0.5)  +  1; 

Vc(l,  N_Spec  *  ( j — 1 )  +  inc_spec)  =  1; 

end; 

end; 
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°/0  Randomly  remove  items  until  solution  is  feasible 
feasible  =  0; 

while  "feasible; 

if  (Vc  *  UPCWV ( : , 2)  <=  Power(i))  &  ... 

(Vc  *  UPCWV ( : ,3)  <=  Cost (i) )  &  ... 

(Vc  *  UPCWV ( : ,4)  <=  Weight (i) )  &  ... 

(Vc  *  UPCWV ( :  ,  5)  <=  Volume (i)); 
feasible  =  1; 

end; 

7,  If  starting  solution  is  not  feasible,  randomly  remove  objects 
if  feasible  ==  0; 
flipped  =  0; 
while  flipped  ==  0 

rand_indx  =  round (  (N_Type  *  N_Spec)  *  rand (1,1)  -  0.5)  +  1 
if  Vc(l ,rand_indx)  ==  1; 

Vc(l ,rand_indx)  =  0; 
flipped  =  1; 

end; 

end; 

end; 

end; 

Initial_Soln  =  Vc; 

7,  Begin  actual  simulated  annealing  heuristic 
7,  Set  simulated  annealing  parameters 
T  =  500;  t  Initial  temperature 
r  =  .05;  l  Initialize  cooling  rate 
T_min  =  .01;  7*  Terminal  temperature 

N_Max  =  5;  l  No.  solutions  explored  at  each  temperature 

while  T  >  T_min; 
n  =  0; 

while  n  <=  N_Max; 

7,  Generate  random  neighbor  of  current  solution 
feasible  =  0; 
loop_count  =  0; 

while  "feasible  &  loop_count  <=  100; 

Vn  =  Neighbor (Vc ,N_Type ,N_Spec) ; 
if  (Vn  *  UPCWV ( : , 2)  <=  Power(i))  &  ... 
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(Vn  *  UPCWV ( : ,3)  <=  Cost(i))  &  ... 

(Vn  *  UPCWV ( : ,4)  <=  Weight (i))  &  ... 

(Vn  *  UPCWV ( : , 5)  <=  Volume (i) ) ; 
feasible  =  1; 

end; 

loop_count  =  loop_count+l ; 

end; 

l  Compare  current  solution  to  neighboring  solution 
Vc_Util  =  Vc  *  UPCWV(:,1); 

Vn_Util  =  Vn  *  UPCWV(:,1); 

l  If  neighbor  is  better,  move  to  it 
if  (Vn_Util  >  Vc_Util)  &  feasible  ==  1 
Vc  =  Vn; 

l  If  neighbor  is  worse  move  to  it  with  a  probability 
else 

if  rand(l , 1)  <  ... 

exp((Vn_Util  -  Vc_Util)/T)  &  feasible  ==  1 
Vc  =  Vn; 

end; 

end; 

n  =  n+1; 

end; 

7,  Decrease  temperature 
T  =  (1-r)  *  T; 

end; 

/  Compute  total  bus  utility 
Utility (l,i)  =  Vc  *  UPCWV(:,1); 

7#  Update  numbers  of  each  payload  type  in  constellation 
for  j  =  1  :  N_Type; 

for  k  =  1  :  N_Spec; 

indx  =  N_Spec  *  (j-1)  +  k; 
if  Vc(l,indx)  ==  1; 

Soln(i,j)  =  k; 

for  1  =  1  :  PD(k,2,j)+l; 

ENum(j ,Nl(i)+l)  =  ENum(j ,Nl(i)+l)  ... 

+  Death_Exp(Death_Coef ,PD(k,2, j) ,1-1) ; 

end; 
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262  °/0  Sum  utility  of  all  satellites,  and  store  resulting  utility  along  with 

263  °/0  specif icaitons 

264  tot_util  =  0; 

265  for  i  =  1  :  N_Sat ; 

266  tot_util  =  tot_util  +  Utility(l ,i) ; 

267  end; 

268  rep_soln( : , : , replicate)  =  Soln; 

269  rep_utility (replicate)  =  tot_util; 

270  end; 

271 

272  l  Find  replication  with  greatest  utility 

273  indx  =  1 ; 

274  for  i  =  1  :  N_replicate 

275  if  rep_utility (i)  >=  rep_utility (indx) ; 

276  indx  =  i; 

277  end; 

278  end ; 

279 

280  /  Stop  clock 

281  t  =  toe; 

282  t 

283 

284  '/.  Print  best  utility  and  specifications 

285  rep_utility (indx) 

286  rep_soln( : , : , indx) 
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Appendix  D.  Random  Neighbor  Function 


y.7.y:/.7.y:/:/.y:/:/.y.y:/;/.y:/;/.7.y;/.7.y:/:/.y:/:/.y.y:/.y.y:/:/.y:/;/.7.y:/.7.y:/:/.y:/:/.y.y:/.y.y:/.y.7.,/:/.7.,/:/.7.,/:/.7.y.y.7.y.y:/.y.y:/:/. 


#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 


7,  March  2006 

7,  This  function,  used  in  the  simulated  annealing  code,  generates  a  random 
7,  neighboring  solution  to  a  given  solution.  A  neighbor  dif feres  by  the 
7,  exclusion  or  inclusion  of  a  single  payload  type/MMD  combination. 


7o7.7o7o7.7o7o7.7o7o7o7o7o7o7o7.7o7o7.7o7o7.7o7o7.7o7o7o7o7o7o7o7o7o7o7o7o7o7.707o7.7o7o7o7o7o7.7o7.7o7o7o707o7.7o7o7.7o7o7.7o7o7o7o7.7«7o7o7o7o7o707o 


function  Neighbor  =  n(Vc ,N_Type ,N_Spec) ; 

Neighbor  =  Vc; 

length  =  size (Neighbor ,2) ; 


7,  Determine  whether  item  will  be  added  (move  =  1)  or  removed  (move  =  0) . 
7,  If  empty,  can  only  add  items 
if  Neighbor  ==  zeros (1 , length) ; 
move  =  1; 

else 

7,  Check  to  see  one  payload  of  each  type  has  been  included.  In  this 
7.  case,  no  more  items  can  be  added 
full  =  1; 

for  i  =  1  :  N_Type 
type_inc  =  0; 
for  j  =  1  :  N_Spec 

indx  =  N_Spec  *  (i-1)  +  j ; 
if  Neighbor (1 , indx)  ==  1; 
type_inc  =  1; 

end; 

end; 

full  =  full  &  type_inc; 

end; 


7,  If  full,  items  can  only  be  removed,  not  added 
if  full  ==  1; 

move  =  0; 
else ; 

7,  If  neither  empty  or  full,  randomly  decide  to  add/remove  item 
move  =  round (rand (1 , 1) ) ; 

end; 

end; 
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42  l  If  adding  an  object 

43  if  move  ==  1; 

44  °/*  Determine  random  type  to  include  that  is  not  already  included 

45  type_inc  =  1 ; 

46  while  type_inc  ==  1; 

47  type  =  round(N_Type  *  rand(l,l)  -  0.5)  +  1; 

48  type_incl  =  0; 

49 

50  for  i  =  1  :  N_Spec; 

51  if  Neighbor (1 ,N_Spec  *  (type-1)  +  i)  ==  1 ; 

52  type_incl  =  1; 

53  end; 

54  end ; 

55  type_inc  =  type_inc  &  type_incl; 

56 

57  if  type_inc  ==  0; 

58  spec  =  round (N_Spec  *  rand (1,1)  -  0.5)  +  1; 

59  Neighbor (1,  N_Spec  *  (type-1)  +  spec)  =  1; 

60  end ; 

61  end; 

62  “/*  If  removing  object 

63  else 

64  °/«  Determine  random  type  to  remove  that  is  not  already  removed 

65  type_inc  =  0; 

66  while  type_inc  ==  0; 

67  type  =  round (N_Type  *  rand (1,1)  -  0.5)  +  1; 

68  type_incl  =  0; 

69 

70  %  Determine  if  type  is  included 

71  for  i  =  1  :  N_Spec; 

72  if  Neighbor (1 ,N_Spec  *  (type-1)  +  i)  ==  1 ; 

73  type_incl  =  1; 

74  spec  =  i; 

75  end ; 

76  end ; 

77  type_inc  =  type_inc  I  type_incl; 

78 

79  if  type_inc  ==  1; 

80  Neighbor (1,  N_Spec  *  (type-1)  +  spec)  =  0; 

81  end; 

82  end ; 

83  end ; 
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18 
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Appendix  E.  5-Norm  Heuristic  Code 


y.7.y:/.7.y:/:/.y:/:/.y.y:/;/.y:/;/.7.y;/.7.y:/:/.y:/:/.y.y:/.y.y:/:/.y:/;/.7.y:/.7.y:/:/.y//:/.y.y:/.y.y:/.y.7.,/:/.7.,/:/.7.,/:/.7.y.y.7.y.y:/.y.y:/:/. 


#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 


7,  March  2006 

7,  This  program  applies  the  5-norm  heuristic  to  the  selection  and 
7,  specification  of  payloads  for  a  multi-satellite  constellation.  The 
7,  heuristic  is  applied  to  each  of  the  satellites  in  succession  seeking  to 
7,  maximize  the  utility  of  each.  The  heuristic  solves  the  multi-choice, 

7,  multidimensional  knapsack  problem  associated  with  each  bus  through  an 
7,  extension  of  the  traditional  profit-to-cost  ratio  heuristic  for  the 
7,  one-dimensional  knapsack  problem. 


7o7.7o7o7.7o7o7.7o7o7.7o7«707o7o707o7.7o7o7.7o7o7.7o7o7.7o7«7.7o7o7o7o7o7o7o7.7o7o7.7o7o7.7o7o7.7o7o7o7o7o7o7o7.7o7o7.707o7.7o7o7.7o7o7.7o7o7o7o7o7o7o 


clear  all; 

7,  Start  clock 
tic ; 

7,  Set  norm  to  be  the  5-norm 
norm  =  5; 

7,  Set  survival  and  utility  function  tuning  parameters 
Death_Coef  =  log(.9); 

Util_Coef  =  log(.5); 

7,  Set  utility  dependence  parameter 
Type_Dep  =  .5; 

7,  Set  number  of  satellites,  payload  types,  and  nonzero  MMD  specifications 
N_Sat  =  7; 

N_Type  =  8; 

N_Spec  =  3; 

7*  Set  time  horizon  of  problem 
Epochs  =  30; 

7,  Set  launch  time  periods 
N1  =  [2  5  7  9  10  11  13  15]  ; 

7,  Input  Payload  Data  —  [Importance,MMD,Power,Cost,Weight, Volume] 

PD(1 , : , 1)  =  [10  3  500  425  450  15]; 


E-l 


42 

PD(2, : 

II 

o 

6 

500 

460 

475 

17]; 

43 

PD (3, : 

,1)  =  [10 

10 

500 

500 

500 

20]  ; 

44 

PD(1 , : 

to 

II 

00 

01 

3 

475 

375 

400 

16]  ; 

45 

PD (2,  : 

,2)  =  [8.5 

6 

475 

405 

415 

18]  ; 

46 

PD (3, : 

to 

II 

00 

01 

10 

475 

430 

420 

19.5] 

47 

PD(1 , : 

,3)  =  [7.5 

3 

425 

410 

430 

10]  ; 

48 

PD (2,  : 

CO 

II 

Ol 

6 

425 

460 

480 

13]  ; 

49 

PD (3, : 

CO 

II 

Ol 

10 

425 

480 

495 

14]; 

50 

PD(1 , : 

,4)  =  [7 

3 

260 

300 

230 

10]  ; 

51 

PD (2,  : 

,4)  =  [7 

6 

260 

350 

280 

13]  ; 

52 

PD (3, : 

,4)  =  [7 

10 

260 

370 

300 

14]; 

53 

PD(1 , : 

,  5)  =  [6 

3 

225 

370 

380 

13]  ; 

54 

PD (2,  : 

,  5)  =  [6 

6 

225 

400 

390 

15.5]  ; 

55 

PD (3, : 

,5)  =  [6 

10 

225 

410 

395 

17.5] 

56 

PD(1, : 

,6)  =  [5.5 

3 

300 

280 

240 

8]; 

57 

PD (2,  : 

,6)  =  [5.5 

6 

300 

320 

290 

9]; 

58 

PD (3, : 

,6)  =  [5.5 

10 

300 

380 

310 

12]  ; 

59 

PD(1 , : 

,7)  =  [5 

3 

275 

150 

280 

7]  ; 

60 

PD (2,  : 

,7)  =  [5 

6 

275 

190 

350 

9.5]  ; 

61 

PD (3, : 

,7)  =  [5 

10 

275 

240 

410 

14]; 

62 

PD(1 , : 

,  8)  =  [3 

3 

175 

270 

225 

4] 

63 

PD (2,  : 

,  8)  =  [3 

6 

175 

310 

260 

5.5]  ; 

64 

PD (3, : 

,  8)  =  [3 

10 

175 

335 

300 

8]; 

65 

66 

7,  Set 

satellite 

resource 

capacities 

67 

for  i 

=  1 :W_Sat ; 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 


Power(i)  =  10500; 

Cost(i)  =  2500; 

Weight (i)  =  2500; 

Volume(i)  =  100; 

end; 

ENum_0rig  =  zeros (N_Type , Epochs) ; 

7.  Set  number  of  payloads  of  each  type  and  specification  initially  in 

#/o  the  constellation 

Init_Cons  =  zeros (N_Type ,N_Spec) ; 

Init.Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0] 

7,  Calculate  expected  number  of  payloads  at  each  time  period 
for  i  =  1: Epochs; 

for  j  =  l:N_Type; 

Num=0 ; 


E-2 


85 

for  k  =  l:N_Spec 

86 

if  i  <=  PD(k,2, j)  +  1; 

87 

Num  =  Num  +  Init_Cons ( j ,k)  *  ... 

88 

Death_Exp(Death_Coef ,PD(k,2, j) ,i-l) ; 

89 

end; 

90 

end; 

91 

ENum(j,i)  =  Num; 

92 

end; 

93 

end; 

94 

95 

7,  Create  vectors  for  total  satellite  utility  and  payload  specifications 

96 

Utility  =  zeros (1 ,N_Sat) ; 

97 

Soln  =  zeros (N_Sat ,N_Type) ; 

98 

99 

for  i=l:N_Sat; 

100 

101 

7.  Create  vectors  to  store  the  utility,  power,  cost,  weight,  and  volumes 

102 

7,  of  each  payload  type/MMD  specification 

103 

utility  =  zeros(l,  N_Type  *  N_Spec) ; 

104 

utility_ratio  =  zeros(l,  N_Type  *  N_Spec) ; 

105 

power  =  zeros (1,  N_Type  *  N_Spec) ; 

106 

cost  =  zeros(l,  N_Type  *  N_Spec) ; 

107 

weight  =  zeros(l,  N_Type  *  N_Spec) ; 

108 

volume  =  zeros(l,  N_Type  *  N_Spec) ; 

109 

110 

for  j  =  1  :  N_Type; 

111 

for  k  =  1  :  N_Spec; 

112 

113 

7,  Load  cost,  weight,  volume,  importance,  and  MMD  of  the  payload 

114 

%  type/MMD  combination 

115 

cost(l,  N_Spec  *  ( j  —  1 )  +  k  )  =  PD(k,4,j); 

116 

weight (1,  N_Spec  *  (j-1)  +  k  )  =  PD(k,5,j); 

117 

volume (1,  N_Spec  *  (j-1)  +  k  )  =  PD(k,6,j); 

118 

psi  =  PD(k, 1 , j ) ; 

119 

MMD  =  PD(k,2, j) ; 

120 

121 

7.  Initialize  total  utility  and  power  consumption  to  zero 

122 

u  =  0; 

123 

ul  =  0; 

124 

p  =  0; 

125 

126 

7,  Compute  power  scaling  factor 

127 

a_l  =  PD(k,3,j)  /  Util_Exp(psi , 1 ,Type_Dep,Util_Coef ,MMD,0) ; 

E-3 

128 


7,  Compute  total  utility  and  power  consumption  of  a  combination 
for  1  =  1  :  MMD  +  1; 


129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 


l  Compute  expected  number  if  combination  is  included 
N  =  max(  ENum(j ,Nl(i)+l)  +  ... 

Death_Exp(Death_Coef ,  MMD,  1-1),  1); 
ul  =  Util_Exp(  psi,N,Type_Dep,Util_Coef , MMD, 1-1  )  ... 

*  Death_Exp(  Death_Coef , MMD, 1-1  ); 
u  =  u  +  ul; 
p  =  p  +  ul  *  a_l; 

end; 

7.  Update  utility  and  power  arrays  with  computed  values 
utility(l,  N_Spec  *  (j-1)  +  k  )  =  u; 
power(l,  N_Spec  *  (j-1)  +  k  )  =  p; 

7*  Compute  5-norm  of  resource  ratios 

size  =  (  (p/Power (i) ) "norm  +  (PD(k,4,j ) /Cost (i) ) "norm  +  ... 

(PD (k, 5, j) /Weight (i) ) "norm  +  ... 

(PD(k,6, j) /Volume (i)) "norm  )" (1/norm) ; 

7,  Update  vector  of  utility-to-aggregated  weight  ratios 
utility_ratio(l ,  N_Spec  *  (j-1)  +  k  )  =  u/size; 

7,  Create  utility/power/cost/weight/volume  master  array 

UPCWV  =  [utility_ratio’ , utility’ , power’ , cost ’, weight 5 , volume’] ; 

end; 

end; 

7,  After  all  values  for  payloads  have  been  computed 
7,  determine  payloads  to  include 

7,  Sort  the  utility-to-aggregrate  weight  ratios 
utility_ratio  =  sort (utility_ratio) ; 

7,  Initialize  vector  of  combinations  included 
Types_Inc  =  zeros (1 ,N_Type) ; 
for  j  =  1  :  N_Type  *  N_Spec 

u  =  utility_ratio(N_Type*N_Spec+l-j) ; 
ind  =  0; 
indx  =  0; 


7,  Determine  index  of  combination  in  unsorted  UPCWV  master  array 


E-4 


171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 
213 


while  ind  ==  0  &  indx  <  N_Type*N_Spec+2 ; 
if  UPCWV ( indx+ 1,1)  ==  u; 
ind  =  1; 

end; 

indx  =  indx  +  1; 

end; 

7.  Determine  payload  type  of  combination 
Type  =  ceil(indx/N_Spec) ; 

7,  Determine  specification  of  combination 
if  mod (indx, N_Spec)  ==  0; 

Spec  =  N_Spec; 

else 

Spec  =  mod ( indx, N_Spec ) ; 

end; 

7,  Calculate  remaining  resources  if  payload  is  included 
if  Types_Inc(l ,Type)  ==  0; 

cl  =  Power(i)  -  UPCWV(indx,  3); 
c2  =  Cost (i)  -  UPCWV (indx,  4); 
c3  =  Weight (i)  -  UPCWV(indx,  5); 
c4  =  Volume(i)  -  UPCWV(indx,  6); 

7.  If  including  combination  is  feasible,  include  it 
if  cl  >=  0  &  c2  >=  0  &  c3  >=  0  &  c4  >=  0; 

Utility(i)  =  Utility(i)  +  UPCWV(indx,  2); 

Power(i)  =  Power(i)  -  UPCWV(indx,  3); 

Cost(i)  =  Cost(i)  -  UPCWV(indx,  4); 

Weight(i)  =  Weight(i)  -  UPCWV(indx,  5); 

Volume(i)  =  Volume(i)  -  UPCWV(indx,  6); 
Types_Inc(Type)  =  PD(Spec,2,Type) ; 

7#  Update  number  of  payloads  in  constellation 
for  1=1:  Types_Inc(Type)+l 

ENum(Type ,N1 (i)+l)  =  ENum(Type,Nl(i)+l)  +  ... 
Death_Exp(Death_Coef ,  Types_Inc (Type) ,  1-1) 

end; 

end; 

end; 

end; 

7,  Update  payload  specifications  for  satellite 


E-5 


214  Soln(i,:)  =  Types_Inc; 

215  end; 

216 

217  °/0  Display  payload  specifications  on  all  satellites 

218  Soln 

219 

220  %  Stop  clock 

221  t  =  toe; 

222  t 

223  tot_util  =  0; 

224 

225  °f  Compute  total  constellation  utility 

226  for  i  =  l:N_Sat; 

227  tot_util  =  tot_util  +  Utility(l , i) ; 

228  end; 

229 

230  7,  Display  total  utility  and  remaining  resources  of  satellites 

231  tot_util 

232  Power  (1,:) 

233  Cost (1 , : ) 

234  Weight (1,:) 

235  Volume (1 , : ) 


Appendix  F.  Weighted  Norm  Heuristic  Code 


y.7.y;/.7.y:/:/.y:/:/.y.y:/;/.y:/;/.7.y;/.7.y:/:/.y:/:/.y.y:/.y.y:/;/.y////.7.y:/.7.y:/:/.y//:/.y.y:/;/.y:/;/.%y;/.7.y//:/.y.y.7.y.y:/.y.y:/.y.y.,/.y. 

#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 

7,  March  2006 

7,  This  program  applies  the  weighted  heuristic  to  the  payload  selection  and 
7,  specification  problem.  The  weighted  norm  heuristic  applies  a  weighted, 

7.  2-norm  to  the  vector  of  resource  requirement  to  capacity  ratios.  The 
7,  terms  of  the  2-norm  are  weighted  by  the  relative  scarcity  of  each 
7,  resource . 

7o7«7o7o7.7«7«7«7»7.7»7o7«7o7o7«7o7o7.7o7o7o7o7.7o7o7o7»7o7.7«7o7.7«7o7.707o7»7»7o7«7»7o7«7o7«7«7o7.7o7o7.707o7o7o7o7.7»7o7o7»7.7o7o7.7«7o7.7o7o7«7»7o 

clear  all; 

7,  Start  clock 
tic ; 

7,  The  heuristic  uses  a  2-norm 
norm  =  2; 

7.  Set  survival  and  utility  decay  function  tuning  parameters 
Death_Coef  =  log(.9);  Util_Coef  =  log(.5); 

7.  Set  dependence  parameter 
Type_Dep  =  .5; 

7,  Set  number  of  satellites,  payload  types,  and  nonzero  MMD  specifications 
N_Sat  =  7;  N_Type  =  8;  N_Spec  =  3; 

7,  Set  time  horizion  of  problem 
Epochs  =  30; 

7,  Set  launch  time  periods 


32 

N1  = 

[2  5  7  9  10 

11  13 

15]  ; 

33 

34 

7,  Input  Payload 

Data  — 

[Import ance, MMD, Power , Cost , Weight , Volume] 

35 

PD(1, 

: , 1)= [10 

3 

500 

425 

450 

15]  ; 

36 

PD  (2, 

: , 1)= [10 

6 

500 

460 

475 

17]; 

37 

PD  (3, 

:  ,  1)  =  [10 

10 

500 

500 

500 

20]  ; 

38 

PD(1 , 

to 

II 

00 

CJ1 

3 

475 

375 

400 

16]  ; 

39 

PD  (2, 

:  ,2)  =  [8.5 

6 

475 

405 

415 

18]  ; 

40 

PD  (3, 

to 

II 

00 

CJ1 

10 

475 

430 

420 

19.5]  ; 

41 

PD(1 , 

:  ,  3)  =  [7 . 5 

3 

425 

410 

430 

10]  ; 

F-l 


42 

PD(2, : 

,3)  =  [7 .5 

6 

425 

460 

480 

13]  ; 

43 

PD (3, : 

,3)  =  [7.5 

10 

425 

480 

495 

14]; 

44 

PD(1 , : 

,4)  =  [7 

3 

260 

300 

230 

10]; 

45 

PD (2,  : 

,4)  =  [7 

6 

260 

350 

280 

13]  ; 

46 

PD (3, : 

s4)  =  [7 

10 

260 

370 

300 

14]; 

47 

PD(1 , : 

,5)  =  [6 

3 

225 

370 

380 

13]  ; 

48 

PD (2,  : 

,5)  =  [6 

6 

225 

400 

390 

15.5]  ; 

49 

PD (3, : 

,5)  =  [6 

10 

225 

410 

395 

17.5]  ; 

50 

PD(1 , : 

,6)  =  [5 . 5 

3 

300 

280 

240 

8]; 

51 

PD (2,  : 

,6)  =  [5.5 

6 

300 

320 

290 

9]; 

52 

PD (3, : 

,6)  =  [5.5 

10 

300 

380 

310 

12]  ; 

53 

PD(1, : 

,7)  =  [5 

3 

275 

150 

280 

7]  ; 

54 

PD (2,  : 

,7)  =  [5 

6 

275 

190 

350 

9.5]  ; 

55 

PD (3, : 

,7)  =  [5 

10 

275 

240 

410 

14]; 

56 

PD(1 , : 

,  8)  =  [3 

3 

175 

270 

225 

4]; 

57 

PD (2,  : 

,  8)  =  [3 

6 

175 

310 

260 

5.5]  ; 

58 

PD (3, : 

,  8)  =  [3 

10 

175 

335 

300 

8]; 

59 

60 

7,  Set 

satellite 

resource 

capacities 

61  Power_Limit  =  10500; 

62  Cost_Limit  =  2500; 

63  Weight_Limit  =  2500; 

64  Volume_Limit  =  100; 

65 

66  ENum_0rig  =  zeros (N_Type , Epochs) ; 

67 

68  7,  Input  number  of  payload  type/MMD  specification  combinations  in 

69  7,  the  constellation  at  time  0 

70  Init_Cons  =  zeros (N_Type ,N_Spec) ; 

71  Init_Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0] ; 

72 

73  7*  Calculate  expected  number  of  payloads  at  each  time  period 

74  for  i  =  1: Epochs; 

75  for  j  =  l:N_Type; 

76  Num=0 ; 

77  for  k  =  l:N_Spec 

78  if  i  <=  PD(k,2,j)  +  1; 

79  Num  =  Num  +  Init_Cons ( j ,k)  *  ... 

80  Death_Exp(Death_Coef ,PD(k,2, j) ,i-l) ; 

81  end; 

82  end ; 

83  ENum(j,i)  =  Num; 

84  end ; 
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85 

end; 

86 

87 

l  Set  satellite  resource  capacities 

88 

for  i  =  l:N_Sat; 

89 

Power (i)  =  Power_Limit ; 

90 

Cost(i)  =  Cost_Limit; 

91 

Weight (i)  =  Weight_Limit ; 

92 

Volume (i)  =  Volume_Limit ; 

93 

end; 

94 

95 

i  Create  vector  to  store  resource  scarcity  values 

96 

PCWV_Sc  =  zeros (N_Sat ,4) ; 

97 

98 

Z  Create  vector  to  store  total  utility  of  each  satellite  and  its  payload  specifications 

99 

Utility  =  zeros (1 ,N_Sat) ; 

100 

Soln  =  zeros (N_Sat ,N_Type) ; 

101 

102 

for  i  =  l:N_Sat; 

103 

104 

°/«  Create  vectors  to  store  the  utility,  power,  cost,  weight,  and  volumes 

105 

utility  =  zeros(l,  N_Type  *  N_Spec) ; 

106 

utility_ratio  =  zeros(l,  N_Type  *  N_Spec) ; 

107 

power  =  zeros (1,  N_Type  *  N_Spec) ; 

108 

cost  =  zeros (1,  N_Type  *  N_Spec) ; 

109 

weight  =  zeros(l,  N_Type  *  N_Spec) ; 

110 

volume  =  zeros(l,  N_Type  *  N_Spec) ; 

111 

112 

for  j  =  1  :  N_Type; 

113 

for  k  =  1  :  N_Spec; 

114 

115 

7,  Load  combination  resouce  requirements,  importance,  and  MMD 

116 

cost(l,  N_Spec  *  (j-1)  +  k  )  =  PD(k,4,j); 

117 

weight(l,  N_Spec  *  (j-1)  +  k  )  =  PD(k,5,j); 

118 

volumed,  N_Spec  *  (j-1)  +  k  )  =  PD(k,6,j); 

119 

psi  =  PD(k, 1 , j ) ; 

120 

MMD  =  PD(k,2, j) ; 

121 

122 

7,  Initialize  total  utility  and  power  consumption  values  to  zero 

123 

u  =  0; 

124 

ul  =  0; 

125 

p  =  0; 

126 

127 

7,  Compute  power  scaling  factor 
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a_l  =  PD(k,3,j)  /  Util_Exp(psi , 1 ,Type_Dep ,Util_Coef ,MMD,0) ; 
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7#  Compute  payload  type/MMD  specification’s  total  utility  and  power  consumption 
for  1  =  1  :  MMD  +  1; 

7.  Compute  expected  number  if  combination  is  included 
N  =  max(  ENum(j ,Nl(i)+l)  +  ... 

Death_Exp(Death_Coef ,  MMD,  1-1),  1); 
ul  =  Util_Exp(  psi,N,Type_Dep,Util_Coef , MMD, 1-1  )  *  ... 

Death_Exp(  Death_Coef ,MMD , 1-1  ); 
u  =  u  +  ul; 
p  =  p  +  ul  *  a_l; 

end; 

7.  Update  utility  and  power  arrays  with  computed  power  and  utility 
utility(l,  N_Spec  *  (j-1)  +  k  )  =  u; 
power(l,  N_Spec  *  (j-1)  +  k  )  =  p; 


end; 

end; 

PCWV_sc  =  zeros (1,4); 

7,  Compute  total  resource  requirements  of  all  combinations  on  a  bus 
for  j  =  1  :  N_Type; 

for  k  =  1  :  N_Spec; 

PCWV_sc (1 , : )  =  PCWV_sc (1 , : )  +  ... 

[power(l,  N_Spec  *  (j-1)  +  k  ),  ... 

PD(k,4,j),  PD(k,5, j) ,  PD(k,6,j)] ; 

end; 

end; 

7,  Compute  resource  scarcities 
PCWV_sc (1 , 1)  =  PCWV_sc (1 , 1)  /  Power (i); 

PCWV_sc (1 , 2)  =  PCWV_sc (1 , 2)  /  Cost(i); 

PCWV_sc(l ,3)  =  PCWV_sc (1 , 3)  /  Weight(i); 

PCWV_sc (1 ,4)  =  PCWV_sc(l ,4)  /  Volume(i); 

7,  Compute  weighted  norm  and  prof it-to-requrements  ratios  of  all  combinations 
for  j  =  1  :  N_Type; 

for  k  =  1  :  N_Spec; 

size  =  (  PCWV_sc  (1 , 1)  *  (p/Power  (i)  ) ''norm  +  ... 

PCWV_sc(l,2)*(PD(k,4, j ) /Cost (i) ) "norm  +  ... 


F-4 


171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 
213 


PCWV_sc(l,3)*(PD(k,5, j ) /Weight (i) ) “norm  +  ... 

PCWV_sc (1,4)* (PD (k , 6 , j ) /Volume (i) ) “norm  ) “ (1/norm) ; 
utility_ratio(l ,  N_Spec  *  (j-1)  +  k  )  =  ... 
utility (1 ,N_Spec  *  (j-1)  +  k  )  /  size; 

%  Create  utility/power/cost/weight/volume  master  array 

UPCWV  =  [utility_ratio’ , utility’ , power’ , cost ’, weight 5 , volume*] ; 

end; 

end; 

°/«  Sort  the  ratios  based  on  the  weighted  norm 
utility_ratio=sort (utility_ratio) ; 

°/0  Initialize  vector  of  included  payload  types 
Types_Inc  =  zeros (1 ,N_Type) ; 

°/0  Iterate  over  all  combinations 
for  j  =  1  :  N_Type  *  N_Spec 

u  =  utility_ratio(N_Type*N_Spec+l-j) ; 
ind  =  0; 
indx  =  0; 


/  Determine  index  of  combination  in  unsorted  UPCWV  master  array 
while  ind  ==  0  &  indx  <  N_Type*N_Spec+2 ; 
if  UPCWV (indx+ 1,1)  ==  u; 
ind  =  1; 

end; 

indx  =  indx  +  1; 

end; 

/  Determine  payload  type  and  MMD  specification 
Type  =  ceil(indx/N_Spec) ; 
if  mod (indx, N_Spec)==0; 

Spec  =  N_Spec; 

else 

Spec  =  mod(indx,N_Spec) ; 

end; 

/  Calculate  remaing  resources  if  combination  is  included 
if  Types_Inc(l ,Type)  ==  0; 

cl  =  Power(i)  -  UPCWV(indx,  3); 
c2  =  Cost (i)  -  UPCWV (indx,  4); 
c3  =  Weight (i)  -  UPCWV(indx,  5); 
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c4  =  Volume(i)  -  UPCWV(indx,  6); 
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234 

235 

236 

237 
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7o  If  including  combination  is  feasible,  include  it 
if  cl  >=  0  &  c2  >=  0  &  c3  >=  0  &  c4  >=  0; 

Utility(i)  =  Utility(i)  +  UPCWV(indx,  2); 

Power(i)  =  Power(i)  -  UPCWV(indx,  3); 

Cost(i)  =  Cost(i)  -  UPCWV(indx,  4); 

Weight(i)  =  Weight(i)  -  UPCWV(indx,  5); 

Volume(i)  =  Volume(i)  -  UPCWV(indx,  6); 

Types_Inc(Type)  =  PD (Spec, 2, Type) ; 

'/.  Update  expected  number  of  payloads  in  constellation  at  each  time  period 
for  1=1:  Types_Inc(Type)+l 

ENum (Type , N1 (i) +1)  =  ENum(Type,Nl(i)+l)  +  ... 

Death_Exp(Death_Coef ,  Types_Inc (Type) ,  1-1); 

end; 

end; 

end; 

end; 

°/*  Store  payload  specifications 
Soln(i,:)  =  Types_Inc; 

end; 

°/0  Display  payload  specifications 
Soln 

°/0  Stop  clock 
t  =  toe; 


°/0  Compute  total  constellation  utility 
tot_util  =  0; 
for  i  =  l:N_Sat; 

tot_util  =  tot_util  +  Utility (1 , i) ; 

end; 

°/0  Display  total  constellation  utility  and  remaining  satellite  resources 

tot_util 

Power (1 , : ) 

Cost  (1 ,  :  ) 

Weight (1 , : ) 

Volume (1 , : ) 


Appendix  G.  Greedy  Heuristic  Code 


y.7.y;/.7.y:/:/.y:/:/.y.y:/;/.y:/:/.7.y;/.7.y:/:/.y:/:/.y.7:/.y.y:/;/.y:/;/.7.y:/.7.y:/:/.y//:/.y.y:/;/.y:/;/.%y;/.7.y//:/.y.y.7.y.y:/.y.y:/.y.7.,/.y. 

#/o  AUTHOR:  Capt  John  Flory 
7,  AFIT/ENS/ G0R-06M 

7,  March  2006 

7,  This  program  applies  a  greedy  heuristic  to  each  satellite  in  the  payload 
7,  selection  and  specification  problem.  All  payload  type/MMD  specifications 
7,  are  sorted  by  their  total  utility.  Payloads  are  included  in  order  of 
7,  decreasing  total  utility.  If  the  includsion  of  a  payload  is  feasible, 

7,  the  payload  is  included. 

7o7«7o7o7.7»7«7.7»7.7»7o7«7o7o7«7o7o7.7o7o7.7»7.7o7»7o7»7o7.7o7o7.7«7o7.707o7»7»7o7«7»7o7o7o7«7o7o7«7o7o7.707o7o7o7o7.7»7o7o7»7.7o7o7.7«7o7.7o7o7«7»7o 

7,  Start  clock 
tic ; 

7,  Set  survival  and  utility  decay  function  tuning  parameters 
Death_Coef  =  log(.9);  Util_Coef  =  log(.5); 

7,  Set  utility  dependence  parameter 
Type_Dep  =  .5; 

7,  Set  number  of  satellites,  payload  types,  and  nonzero  MMD  specifications 
N_Sat  =  7;  N_Type  =  8;  N_Spec  =  3; 

7,  Set  problem  time  horizion 
Epochs  =  30; 

7*  Input  launch  time  periods 


28 

N1  = 

[2  5  T  9  10 

11  13 

15]  ; 

29 

30 

7,  Input  Payload 

Data  — 

[Importance ,  MMD ,  Power ,  Cost ,  Weight ,  Volume] 

31 

PD(1, 

II 

o 

3 

500 

425 

450 

15]  ; 

32 

PD  (2, 

II 

o 

6 

500 

460 

475 

17]; 

33 

PD  (3, 

II 

o 

10 

500 

500 

500 

20]  ; 

34 

PD(1 , 

to 

II 

00 

U1 

3 

475 

375 

400 

16]  ; 

35 

PD  (2, 

to 

II 

00 

CJ1 

6 

475 

405 

415 

18]  ; 

36 

PD  (3, 

to 

II 

00 

01 

10 

475 

430 

420 

19.5]  ; 

37 

PD(1 , 

:  ,  3)  =  [7 . 5 

3 

425 

410 

430 

10]  ; 

38 

PD  (2, 

CO 

II 

Ol 

6 

425 

460 

480 

13]  ; 

39 

PD  (3, 

: ,  3)  =  [7 . 5 

10 

425 

480 

495 

14]; 

40 

PD(1, 

II 

3 

260 

300 

230 

10]; 

41 

PD  (2, 

n 

6 

260 

350 

280 

13]  ; 
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42 

PD(3, : 

,4)  = 

[7 

10 

260 

370 

300 

14]; 

43 

PD(1 , : 

s5)  = 

[6 

3 

225 

370 

380 

13]  ; 

44 

PD (2,  : 

s5)  = 

[6 

6 

225 

400 

390 

15.5]  ; 

45 

PD (3, : 

s5)  = 

[6 

10 

225 

410 

395 

17.5]  ; 

46 

PD(1, : 

,6)  = 

[5 

5 

3 

300 

280 

240 

8]; 

47 

PD (2,  : 

,6)  = 

[5 

5 

6 

300 

320 

290 

9]; 

48 

PD (3, : 

s6)  = 

[5 

5 

10 

300 

380 

310 

12]  ; 

49 

PD(1 , : 

,7)  = 

[5 

3 

275 

150 

280 

7]; 

50 

PD (2,  : 

,7)  = 

[5 

6 

275 

190 

350 

9.5]  ; 

51 

PD (3, : 

,7)  = 

[5 

10 

275 

240 

410 

14]; 

52 

PD(1 , : 

,8)  = 

[3 

3 

175 

270 

225 

4]; 

53 

PD (2,  : 

,8)  = 

[3 

6 

175 

310 

260 

5.5]  ; 

54 

PD (3, : 

,8)  = 

[3 

10 

175 

335 

300 

8]; 

55 

56 

7.  Set 

satellite  bus 

resource 

capacities 

57 

for  i 

=  1: 

N_Sat ; 

58  Power(i)  =  10500; 

59  Cost(i)  =  2500; 

60  Weight (i)  =  2500; 

61  Volume(i)  =  100; 

62  end ; 

63 

64  ENum_0rig  =  zeros (N_Type , Epochs) ; 

65 

66  7,  Set  numbers  of  each  payload  type/MMD  specification  in  the  constellation 

67  7,  at  time  0 

68  Init_Cons  =  zeros (N_Type ,N_Spec) ; 

69  Init_Cons  =  [0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0;  0  0  0] ; 

70 

71  7o  Calculate  expected  number  of  payloads  at  each  time  period 

72  for  i  =  1: Epochs; 

73  for  j  =  1 : N_Type ; 

74  Num=0 ; 

75  for  k  =  l:N_Spec 

76  if  i  <=  PD(k,2,j)  +  1; 

77  Num  =  Num  +  Init_Cons ( j ,k)  *  ... 

78  Death_Exp(Death_Coef ,PD(k,2, j) ,i-l) ; 

79  end ; 

80  end ; 

81  ENum(j,i)  =  Num; 

82  end ; 

83  end ; 

84 
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85 

l  Create  vectors  to  store  utility  and  payload  specifications 

86 

Utility  =  zeros (1 ,N_Sat) ; 

87 

Soln  =  zeros (N_Sat ,N_Type) ; 

88 

89 

for  i=l:N_Sat; 

90 

91 

7.  Create  vectors  to  store  the  utility,  power,  cost,  weight,  and  volume 

92 

7.  requirements 

93 

utility  =  zeros(l,  N_Type  *  N_Spec) ; 

94 

power  =  zeros (1,  N_Type  *  N_Spec) ; 

95 

cost  =  zeros(l,  N_Type  *  N_Spec) ; 

96 

weight  =  zeros (1,  N_Type  *  N_Spec) ; 

97 

volume  =  zeros (1,  N_Type  *  N_Spec) ; 

98 

99 

for  j  =  1  :  N_Type; 

100 

for  k  =  1  :  N_Spec; 

101 

7,  Load  combination’s  cost,  weight,  volume,  importance,  and  MMD 

102 

cost(l,  N_Spec  *  (j-1)  +  k  )  =  PD(k,4,j); 

103 

weight(l,  N_Spec  *  (j-1)  +  k  )  =  PD(k,5,j); 

104 

volumed,  N_Spec  *  (j-1)  +  k  )  =  PD(k,6,j); 

105 

psi  =  PD(k, 1 , j ) ; 

106 

MMD  =  PD(k,2, j) ; 

107 

108 

7.  Initialize  total  utility  and  power  consumption  summation 

109 

7,  values  to  zero 

110 

u  =  0; 

111 

ul  =  0; 

112 

p  =  0; 

113 

114 

7,  Compute  power  scaling  factor 

115 

a_l  =  PD(k,3, j )  /  Util_Exp(psi , 1 ,Type_Dep,Util_Coef ,MMD,0) ; 

116 

117 

7.  Compute  total  utility  and  power  consumption  of  a  payload 

118 

7.  type/MMD  specification  combination 

119 

for  1  =  1  :  MMD  +  1; 

120 

7#  Compute  expected  number  of  payloads  if  combination  is 

121 

7.  included 

122 

N  =  max(  ENum( j ,Nl(i)+l)  +  ... 

123 

Death_Exp(Death_Coef ,  MMD,  1-1),  1); 

124 

ul  =  Util_Exp(  psi ,N,Type_Dep,Util_Coef , MMD, 1-1  )  *  ... 

125 

Death_Exp(  Death_Coef ,MMD , 1-1  ); 

126 

u  =  u  +  ul; 

127 

p  =  p  +  ul  *  a_l; 
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end; 

7,  Update  utility  and  power  vectors  with  computed  values 
utility(l,  N_Spec  *  (j-1)  +  k  )  =  u; 
power (1,  N_Spec  *  (j-1)  +  k  )  =  p; 

%  Create  utility/power/cost/weight/volume  master  array 
UPCWV  =  [utility’ , power’ , cost ’ , weight’ , volume ’] ; 

end; 

end; 

7,  Sort  total  utility  values 
utility  =  sort (utility) ; 

7,  Initialize  vector  of  included  payload  types 
Types_Inc  =  zeros (1 ,N_Type) ; 

7,  Iterate  over  all  combinations 
for  j  =  1  :  N_Type  *  N_Spec 

u  =  utility (N_Type*N_Spec+l-j ) ; 
ind  =  0; 
indx  =  0; 


7,  Determine  index  of  utility  in  unsorted  UPCWV  master  array 
while  ind  ==  0  &  indx  <  N_Type*N_Spec+2 ; 
if  UPCWV (indx+ 1,1)  ==  u; 
ind  =  1; 

end; 

indx  =  indx  +  1; 

end; 

7,  Determine  payload  type 
Type  =  ceil(indx/N_Spec) ; 

7,  Determine  payload  specification 
if  mod (indx, N_Spec)==0; 

Spec  =  N_Spec; 

else 

Spec  =  mod(indx,N_Spec) ; 

end; 

7,  Calculate  remaining  satellite  resources  if  payload  is  included 
if  Types_Inc(l ,Type)  ==  0; 

cl  =  Power(i)  -  UPCWV(indx,  2); 
c2  =  Cost (i)  -  UPCWV (indx,  3); 
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c3  =  Weight (i)  -  UPCWV(indx,  4); 
c4  =  Volume (i)  -  UPCWV(indx,  5); 

%  If  payload  can  be  included,  include  it 
if  cl  >=  0  &  c2  >=  0  &  c3  >=  0  &  c4  >=  0; 

Utility(i)  =  Utility(i)  +  u; 

Power(i)  =  Power(i)  -  UPCWV(indx,  2); 

Cost(i)  =  Cost(i)  -  UPCWV(indx,  3); 

Weight(i)  =  Weight(i)  -  UPCWV(indx,  4); 

Volume(i)  =  Volume(i)  -  UPCWV(indx,  5); 

Types_Inc(Type)  =  PD (Spec, 2, Type) ; 

/  Update  expected  number  of  payloads  in  constellation  at 

%  each  time  period 

for  1=1:  Types_Inc(Type)+l 

ENum (Type , N1 (i) +1)  =  ENum(Type,Nl(i)+l)  +  ... 

Death_Exp(Death_Coef ,  Types_Inc (Type) ,  1-1); 

end; 

end; 

end; 

end; 

°/0  Store  payload  specifications 
Soln(i,:)  =  Types_Inc; 

end; 

#/o  Stop  clock 
t  =  toe;  t 

7.  Display  payload  specifications 
Soln 

7,  Compute  total  constellation  utility 
tot_util  =  0; 
for  i  =  l:N_Sat; 

tot_util  =  tot_util  +  Utility (1 , i) ; 

end; 

7,  Display  total  utility  and  remaining  resources  on  all  satellites 

tot_util 

Power (1 , : ) 

Cost  (1 ,  :  ) 

Weight (1 , : ) 

Volume (1 , : ) 
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Appendix  H.  Payload  Survival  Function 


%  AUTHOR:  Capt  John  Flory 
7.  AFIT/ENS/ G0R-06M 


7,  March  2006 

7,  This  function  is  called  by  the  programs  to  give  the  exponential 
7,  survival  distribution  of  satellite  payloads. 
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7,  Function  uses  the  following  parameters: 

7,  alpha  -  Tuning  parameter  >  0 

7,  MMD  -  Mean  mission  duration  specification 

7,  n  -  Time  period 

function  Death_Exp  =  F (alpha, MMD ,n) ; 

Death_Exp  =  exp (-abs (alpha) /MMD*n) ; 
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Appendix  I.  Payload  Utility  Decay  Function 
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7*  This  function  is  called  by  the  programs  to  give  the  exponential 
#/o  decay  function  of  payload  utility. 
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7,  Function  uses  the  following  parameters: 

7*  psi  -  Payload  importance 

7,  N  -  Number  of  like-type  functional  payloads 

7,  Type_Dep  -  Value  of  dependence  parameter  (gamma) 

7,  alpha  -  Tuning  parameter 

7,  MMD  -  Payload  mean  mission  duration 

7,  n  -  Time  period 


function  Util_Exp  =  u (psi ,N,Type_Dep, alpha, MMD, n) ; 


7,  Positive  utility  dependence 
if  Type_Dep  >  0; 

Util_Exp  =  psi/N~ (Type_Dep) *exp (-abs (alpha) /MMD*n) ; 

7,  No  utility  dependence 
elseif  Type_Dep  ==  0; 

Util_Exp  =  psi*exp (-abs (alpha) /MMD*n) ; 


7,  Logarithmic  utility  dependence 
elseif  Type_Dep  ==  -1; 

Util_Exp  =  psi/log(max(N, 1) ) *exp (-abs (alpha) /MMD*n) ; 


end; 
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